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Abstract 

Several 7-ray binaries have been recently detected by the High-Energy Stere- 
oscopy Array (H.E.S.S.) and the Major Atmospheric Imaging Cerenkov (MAGIC) 
telescope. In at least two cases, their nature is unknown. In this paper we aim to 
provide the details of a theoretical model of close 7-ray binaries containing a young 
energetic pulsar as compact object, earlier presented in recent Letters. This model 
includes a detailed account of the system geometry, the angular dependence of pro- 
cesses such as Klein-Nishina inverse Compton and 77 absorption in the anisotropic 
radiation field of the massive star, and a Monte Carlo simulation of leptonic cas- 
cading. We present and derive the used formulae and give all details about their 
numerical implementation, particularly, on the computation of cascades. In this 
model, emphasis is put in the processes occurring in the pulsar wind zone of the 
binary, since, as we show, opacities in this region can be already important for close 
systems. We provide a detailed study on all relevant opacities and geometrical de- 
pendencies along the orbit of binaries, exemplifying with the case of LS 5039. This 
is used to understand the formation of the very high-energy lightcurve and phase 
dependent spectrum. For the particular case of LS 5039, we uncover an interesting 
behavior of the magnitude representing the shock position in the direction to the 
observer along the orbit, and analyze its impact in the predictions. We show that 
in the case of LS 5039, the H.E.S.S. phenomenology is matched by the presented 
model, and explore the reasons why this happens while discussing future ways of 
testing the model. 
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1 Introduction 



Very recently, a few massive binaries have been identified as variable very- 
high-energy (VHE) 7-ray sources. They are PSR B1259-63 (Aharonian et al. 
2005a), LS 5039 (Aharonian ct al. 2005b, 2006), LS I +61 303 (Albert et al. 
2006, 2008a,b), and Cyg X-1 (Albert ct al. 2007). The nature of only two of 
these binaries is considered known: PSR B 1259-63 is formed with a pulsar 
whereas Cyg X-1 is formed with a black hole compact object. The nature of the 
two remaining systems is under discussion. The high-energy phenomenology 
of Cyg X-1 is different from that of the others. It has been detected just once 
in a flare state for which a duty cycle is yet unknown. The three other sources, 
instead, present a behavior that is fully correlated with the orbital period. The 
latter varies from about 4 days in the case of LS 5039 to several years in the 
case of PSR B1259-63: this span of orbital periodicities introduces its own 
complications in analyzing the similarities among the three systems. 

LS I -|-61 303 shares with LS 5039 the quality of being the only two known 
microquasars/7-ray binaries that are spatially coincident with sources above 
100 MeV listed in the Third Energetic Gamma-Ray Experiment (EGRET) 
catalog (Hartman et al. 1999). These sources both show low X-ray emission 
and variability, and no signs of emission lines or disk accretion. For LS 1 +61 
303, extended, apparently precessing, radio emitting structiues at angular ex- 
tensions of 0.01-0.05 arcsec have been reported by Massi et al. (2001, 2004); 
this discovery has earlier supported its microquasar interpretation. But the 
uncertainty as to what kind of compact object, a black hole or a neutron star, 
is part of the system (e.g., Casares et al. 2005a), seems settled for many after 
the results presented by Dhawan et al. (2006). These authors have presented 
observations from a July 2006 VLBl campaign in which rapid changes are 
seen in the orientation of what seems to be a cometary tail at periastron. 
This tail is consistent with it being the result of a pulsar wind. Indeed, no 
large features or high-velocity flows were noted on any of the observing days, 
which implies at least its non-permanent nature. The changes within 3 hours 
were found to be insignificant, so the velocity can not be much over 0.05c. 
Still, discussion is on-going (e.g. see Romero et al. 2007, Zdziarski et al. 2008). 
New campaigns with similar radio resolution, as well as new observations in 
the 7-ray domain have been obtained since the Dhawan's et al. original results 
(Albert et al. 2008b). A key aspect in these high-angular-resolution campaigns 
is the observed maintenance in time of the morphology of the radio emission 
of the system: the changing morphology of the radio emission along the orbit 
would require a highly unstable jet, which details are not expected to be re- 
produced orbit after orbit as indicated by current results (Albert et al. 2008b). 
The absence of accretion signatures in X-rays in Chandra and XMM-Newton 
observations (as reported by Sidoli et al. 2006, Chernyakova et al. 2006, and 
Paredes et al. 2007) is another relevant aspect of the discussion about the 
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compact object companion. Finally, it is interesting to note that neutrino de- 
tection or non-dctcction with ICECUBE will shed light on the nature of the 
7-ray emission irrespective of the system composition (e.g., Aharonian et al. 
2006b, Torres and Halzen 2007). 

For LS 5039, a periodicity in the 7-ray flux, consistent with the orbital timescale 
as determined by Casarcs ct al. (2005b), was found with amazing precision 
(Aharonian et al. 2006). Short timescale variability displayed on top of this 
periodic behavior, both in flux and spectrum, was also reported. It was found 
that the parameters of power-law fits to the 7-ray data obtained in 0.1 phase 
binning already displayed significant variability. Current H.E.S.S. observations 
of LS 5039 (~ 70 hours distributed over many orbital cycles, Aharonian ct al. 
2006) constitute one of the most detailed datasets of high-energy astrophysics. 
Similarly to LS I -|-61 303, the discovery of a jet-hke radio structure in LS 5039 
and the fact of it being the only radio/X-ray source co-localized with a mildly 
variable (Torres ct al. 2001a,b) EGRET detection, prompted a microquasar 
interpretation (advanced already by Paredes et al. 2000). However, the cur- 
rent mentioned findings at radio and VHE 7-rays in the cases of LS I -|-61 303 
(Dhawan et al. 2006, Albert et al. 2006, 2008b) or PSR B1259-63 (Aharonian 
et al. 2005), gave the perspective that all three systems are different real- 
izations of the same scenario: a pulsar-massive star binary. Dubus (2006a,b) 
has studied these similarities. He provided simulations of the extended radio 
emission of LS 5039 showing that the features found in high resolution radio 
observations could also be interpreted as the result of a pulsar wind. Recently, 
Ribo et al. (2008) provided VLBA radio observations of LS 5039 with mor- 
phological and astrometric information at milliarcsecond scales. They showed 
that a microquasar scenario cannot easily explain the observed changes in mor- 
phology. All these results, together with the assessment of the low X-ray state 
(Martocchia et al. 2005) made the pulsar hypothesis tenable, and the possi- 
bility of explaining the H.E.S.S. phenomenology in such a case, an interesting 
working hypothesis. 

High energy emission from pulsar binaries has been subject of study for a 
long time (just to quote a non-exhaustive list of references note the works 
of Maraschi and Treves 1981; Protheroe and Stanev 1987, Arons and Tavani 
1993, 1994; Moskalenko et al. 1993; Bednarek 1997, Kirk et al. 1999, Ball and 
Kirk 2000, Romero et al. 2001, Anchordoqui et al. 2003, and others already 
cited above). LS 5039 has been recently subject of intense theoretical studies 
(e.g., Bednarek 2006, 2007; which we comment on in more detail below, Bosch- 
Ramon et al. 2005; Bottcher 2007; Bottcher and Dermer 2005; Dermer and 
Bottcher 2006; Dubus 2006a,b; Paredes et al. 2006; Khangulyan et al. 2007; 
Dubus et al. 2007). 

In the penultimate paper mentioned in the list above, Khangulyan et al. 
(2007), and contrary to the assumption here, authors assumed a jet structure 
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perpendicular to the orbital plane of the system. The energy spectrum and 
lightcurves were computed, accounting for the acceleration efficiency, the loca- 
tion of the accelerator along the jet, the speed of the emitting flow, the inclina- 
tion angle of the system, as well as specific features related to anisotropic in- 
verse Compton (IC) scattering and pair production. Different magnetic fields, 
affecting Synchrotron emission, and the losses they produced, were also tested 
given a large model parameter space. Authors found a good agreement be- 
tween H.E.S.S. data for some of their models. 

In the last of these papers, Dubus et al. (2007) computed the phase depen- 
dent lightcurve and spectra expected from inverse Compton interactions from 
electrons injected close to the compact object, assumed as a likely rotation- 
powered pulsar. Since the angle at which an observer sees the binary and 
propagating electrons changes with the orbit (see below) , a phase dependence 
of the spectrum is expected, and anisotropic inverse Compton is needed to 
compute it. In general, they found that the lightcurve is a good fit to the ob- 
servations, except at the phases of maximum attenuation where pair cascade 
emission plays a role. Dubus et al. (2007) do not consider cascading in their 
models, as we do here. Without cascading, zero flux is expected at a broad 
phase around periastron, which is not found. This lack of cascading in their 
model also affects the spectra, which are not reproduced well, particularly 
at the superior conjunction broad phases of the orbit. They mentioned that 
both, cascading and/or a change in the slope of the power-law injection for 
the interacting electron distribution could be needed to explain the spectrum 
in these phases, what we explore in detail in this work. 

In order to compute inverse Compton emission from LS 5039, we use, as in 
previous works, leptons interacting with the star photon field. Geometry is 
described there with different levels of detail, what inffuence the results. In 
general, cascading processes were not taken into account, and the goodness of 
fitting the H.E.S.S. data is arguable in most cases, both for the hghtcurve and 
spectrum. 

In none of the papers mentioned above, the theoretical predictions for the 
short timescale spectral variability found by H.E.S.S. in 0.1 phase binning 
was shown and compared with data. We discuss these results from our model 
below. 

In recent Letters (Sierpowska-Bartosik and Torres 2007, 2008) under the as- 
sumption that LS 5039 is composed by a pulsar rotating around an 06.5V 
star in the ~ 3.9 day orbit, we presented the results of a leptonic (for a 
generic hadronic model see Romero et al. 2003) theoretical modeling for the 
high-energy phenomenology observed by H.E.S.S. These works studied the 
lightcurve, the spectral orbital variability in both broad orbital phases and in 
shorter (0.1 phase binning) timescales and have found a complete agreement 
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between H.E.S.S. observations and our predictions. We have also analyzed 
how this model could be tested by Gamma-ray Large Area Space Telescope 
(GLAST), and how much time would be needed for this satellite in order to 
rule the model out in case theory significantly departs from reality. But many 
details of implementation which are not only useful for the case of LS 5039 but 
for all others close massive 7-ray binaries, as well as many interesting results 
concerning the binary geometry, wind termination, opacities to different pro- 
cesses along the orbit of the system, and further testing at the highest energy 
7-ray domain were left without discussion in our previous works. Here, we 
provide these details, together with benchmark cases that are useful to under- 
stand the formation of the very high-energy lightcurve and phase dependent 
spectra. 

The rest of this paper is organized as follows. Next Section introduces the 
model concept and its main properties. It provides a discussion of geometry, 
wind termination, and opacities along the orbit of the system (we focus on 
LS 5039). An accompanying Appendix provides mathematical derivations of 
the formulae used and useful intermediate results that are key for the model, 
but too cumbersome to include them as part of the main text. It also deals 
with numerical implementation, and describes in detail the Monte Carlo sim- 
ulation of the cascading processes. The results follow: Section 3 deals with a 
mono-energetic interacting particle population, and Section 4, with power-law 
primary distributions. Comparison with H.E.S.S. results is made in these Sec- 
tions and details about additional tests are given. Final concluding remarks 
are provided at the end. 



2 Description of the model and its implementation 

Under the assumption that the pulsar in the binary is energetic enough to 
prevent matter from the massive companion from accreting, a termination 
shock is created in the interaction region of the pulsar and donor star winds. 
This is represented in Fig. [T] We focus on the specific case of the binary 
LS 5039, which we use as a testbed all along this paper. The volume of the 
system is separated by the termination shock, which structure depends on 
features of the colliding winds: it may be infiuenced by the anisotropy of the 
winds themselves, the motion of the pulsar along the orbit, turbulences in the 
shock fiow, etc. For simplicity it is assumed here that the winds are radial and 
spherically symmetric, and that the termination of the pulsar wind is an axial 
symmetric structure with negligible thickness. In this general picture there 
are three regions of different properties in the binary: the pulsar wind zone 
(PWZ), the shock (SR), and the massive star wind zone (MSWZ). 

The energy content in the interaction population of particles is assumed as a 
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Fig. 1. Left: Sketch of a close binary, such as the LS 5039 system. P stands for the 
orbital position of periastron, A for apastron, INFC for inferior conjunction, and 
SUPC for superior conjunction. The orbital phase, (p, and the angle to the observer, 
o^obsj are marked. The orbital plane is inclined with respect to the direction of the 
observer (this angle is not marked) . The termination shock created in the interaction 
of the pulsar and the massive star winds is also marked for two opposite phases 
-periastron and apastron- together with the direction to the observer for both 
phases. Right: The physical scenario for high-energy photon production in the PWZ 
terminated by the shock. e~^e~ are injected by the pulsar or by a close-to-the-pulsar 
shock and travel towards the observer, producing Inverse Compton photons, 7, via 
up-scattering thermal photons from the massive star, e. 7-photons can initiate IC 
cascade due to absorption in the same thermal field. The cascade is developing up 
to the termination shock. Electron reaching the shock are trapped there in the local 
magnetic field, while photons propagate further and escape from the binary or are 
absorbed in massive star radiation field. The cascade develops radially following the 
initial injection direction given by aobs- 

fraction of pulsar spin-down power Lgd- In case of young energetic pulsars, this 
power is typically ~ 10^^ — lO^^ergs"^. Propagating pairs up-scatter thermal 
photons from the massive star due to inverse Compton process. For close 
binaries, the radiation field of hot massive stars (type O, Be or WR, having 
typical surface temperatures in the range Tg ~ 10^—10^ K and linear dimension 
Rg ~ IO-R0) dominates along the whole orbit over other possible fields (e.g., 
the magnetic field or the thermal field of the neutron star). This thermal 
radiation field is anisotropic, particularly for e~^e~ injected close to the pulsar 
(the radiation source is misplaced with respect to the electron injection place). 

The high-energy photons produced by pairs can initiate cascades due to sub- 
sequent pair production in absorption (77) process with the same radiation 
field (as sketched in Fig. [1]). We assume that these cascades develop along 
the primary injection direction, i.e., in a one- dimensional way, which is cer- 
tainly justified based on the relativistic velocity of the interacting electrons. 
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This process is followed up to the termination shock unless leptons lose their 
energy before reaching it. Those leptons which propagate to the shock region 
are trapped there by its magnetic field. Radiation from them is isotropised. 
The photons produced in cascades which reach the shock can get through it 
and finally escape from the binary or be absorbed in the radiation field close 
to the massive star, some may even reach the stellar surface. 

In the shock region leptons move along it with velocity ~ c/3. They could 
be re-accelerated and produce radiation via synchrotron (local magnetic field 
from the pulsar side) or inverse Compton scattering (ICS, thermal radiation 
field from the massive star) processes. However, as they are isotropised in 
the local magnetic field, photons are produced in different directions and their 
directionality towards the observer is lost. It was already shown by Sierpowska 
and Bednarek (2005) that in compact binary systems (as an example, the 
parameters of Cyg X-3 were taken by these authors) the radiation processes 
in the shock region do not dominate: the energy carried by e^e~ reaching the 
shock is a small fraction of total injected power. Furthermore, we will show 
that for the parameters relevant to the LS 5039 scenario, the PWZ is relatively 
large with respect to the whole volume of the system for the significant range 
of the binary orbit. 



2.1 Hydrodynamic balance 



Assuming that both, the pulsar and the massive star winds are spherically 
symmetric, and based on the hydrodynamic equilibrium of the fiows, the ge- 
ometry of the termination shock is described by parameter rj — MiVi/MoVo, 
where AfjVi and MoVq are the loss mass rates and velocities of the two winds 
(Girard and Wilson, 1987). The shock will be symmetric with respect to the 
line joining two stars, with a shock front at a distance from the one of the 
stars: 

rs = D—^. (1) 



The surface of the shock front can be approximated then by a cone-like 
structure with opening angle given by = 2.1 ^1 — ^-^^ f]^^^, where rj — 
min{ri, r]^^) takes the smaller value between the two magnitudes quoted. This 
last expression was achieved under the assumption of non-relativistic winds 
in the simulations of the termination shock structure in Girard and Wilson 
(1987), albeit it is also in agreement with the relativistic winds case (e.g., 
Eichler and Usov, 1993; Bogovalov et al., 2007). If one of the stars is a pulsar 
of a spin down luminosity Lg^ and the power of the massive star is M^V^, the 
parameter rj can be calculated from the formula rj = ^^^"y ^ (e.g.. Ball and 
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Kirk 2000). Note that for 77 < 1, the star wind dominates over the pulsar's and 
the termination shock wraps around it. Note also that for r/ = 1, the shock 
is at equal distance, (i/2, between the stars. In the case of LS 5039, for the 
assumed (nominal) spin-down luminosity Lgd discussed below, the value of 77 
is between 0.5 (periastron) and 0.3 (apastron). 

The massive star in the LS 5039 binary system is of O type, which wind 
is radiation driven. The velocity of the wind at a certain distance can be 

described by classical velocity law K(r) = Vq + (Voo — Vq) {l — , where 
Voo is the wind velocity at infinity, Rg the hydrostatic radius of the star, Vq 
is the velocity close to the stellar surface, and (5 — 0.8 — 1.5 (e.g., Cassinelli 
1979, Lamers and Cassinelli 1999) and we assume (3 = 1.5. As can be seen 
by plotting the velocity law, the influence of this parameter is minor. Typical 
wind velocities for 0/Be type stars are V^o ~ (1 — 3) x 10"^kms"^. For LS 5039 
we have V^c = 2.4 x 10^ kms"^ and = 4 kms"^ (Casares et al. 2005). The 
assumed value for the star mass- loss rate is 10~^ yr~^, while the typical 
values for 0/Be stars are 10~^ — 10~^ Mq yr~^. 

2.2 The PWZ and the interacting lepton population 

The magnetization parameter a = B"^ / A'lr'-fnmc'^ , (e.g., Langdon et al. 1988) 
is defined as the ratio of Poynting flux to relativistic particle energy flux. 
The magnetic field in a is that of the upstream shock propagating with bulk 
Lorentz factor 7 and n being the relativistic particle density. The processes es- 
tablishing the effective change of a along the PWZ are the central issue in the 
discussion of dissipation mechanisms in relativistic plasma fiows, exemplified 
with the Crab pulsar wind (e.g.. Kennel and Coroniti, 1984), where variations 
seem notable. The Crab wind is originally Foynting-dominated {a ~ 10^ close 
to the neutron star, e.g. Arons 1979); but it is kinetic-dominated near the 
termination shock (cr ~ 10~^, e.g.. Kennel & Coroniti 1984). The change in 
a is produced as a result of dissipative plasma processes in the FWZ, which 
is characterized by high bulk Lorentz factor (e.g., Melatos 1998). Dissipa- 
tion (plasma processes engaged in the conversion of electromagnetic towards 
particle kinetic energy) in Poynting flux dominated plasma flows can be in 
the form of stochastic/non-stochastic and adiabatic/non-adiabatic processes, 
thermal heating/non- thermal particle generation, and isotropic adiabatic ex- 
pansion/directed bulk acceleration of the plasma flow (Jaroschek et al. 2008). 
The microphysical details are decisive when looking for the type of particle 
energization and their spectra. 

The existence of wisps in the inner structure of the Crab nebula has been 
discussed by, e.g., Lou (1998). The interesting fact of some of them being 
close to the pulsar, apparently well inside the PWZ, was interpreted as being 
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produced by slightly inhomogeneous wind streams, demonstrating that reverse 
fast MHD shocks at various spin latitudes can appear quasi-stationary in space 
when their propagation speeds relative to the pulsar wind are comparable to 
the relativistic outflow. The possibihty of an inhomogeneous wind stream is 
not implausible. Successive radio pulses from a pulsar indeed vary in their 
shapes, eventhough the average pulse is stable. A slower wind stream will be 
eventually caught by a faster one to trigger forward and reverse fast MHD 
shocks inside the PWZ. In this zone, charged particles can be further acceler- 
ated by these turbulences and magnetic reconnection. Lou proposed that an 
isotropised power-law like energy distribution of the electrons thus produced 
help to understand the properties of the changing and brilliant inner nebula. 

A mono-energetic assumption for the distribution of leptons in the PWZ can 
be considered as a first approach to the problem. On one hand, the magne- 
tization parameter may be a function of angle, and although must be very 
small in the equatorial part of the wind (e.g.. Kirk 2006), simulations do not 
favor an angle- independent low value (Komissarov & Lyubarsky 2004). On the 
other hand, Contopoulos & Kazanas (2002) already showed that the Lorentz 
factor of the outfiowing plasma could increase linearly with distance from 
the light cylinder (implying that a decreases inversely proportional to the 
distance). Contopoulos & Kazanas (2002) mentioned that this specific radial 
dependence of the pulsar winds Lorentz factor is expected to have additional 
observational consequences: e.g., Bogovalov & Aharonian (2000) computed 
the Comptonization of soft photons to TeV energies in the Crab through their 
interaction with the expanding MHD wind, while Tavani & Arons (1997) and 
Ball & Kirk (2000) computed the corresponding radiation expected by the 
radio-pulsar Be star binary system PSR B1259-63 through the interaction of 
the relativistic wind with the photon field of the companion in much the same 
way we do here. Indeed, the details in these predictions would be modified, 
as shown by Sierpowska & Bednarek (2004, 2005) should the linear accelera- 
tion model be adopted. Hibschman and Arons (2001) discussed the creation of 
electron-positron cascades in the context of pulsar polar cap acceleration mod- 
els. They computed the spectrum of pairs that would be produced outflowing 
the magnetosphere. They found that the pair spectra should be described by 
a power-law. 

One possibility for the dissipative conversion is established by magnetic re- 
connection processes between anti-parallel magnetic stripes during outwards 
propagation in the PWZ (Lyubarsky and Kirk, 2001; Kirk and Skjaeraasen, 
2003). Kirk (2004) considered acceleration in relativistic current sheets (large 
magnetization parameter, with Alfven speed va — C\f[p I (o'+l) close to c). Re- 
cently, Jaroschek et al. (2008, and see references therein for related work) ad- 
dressed the problem of interacting relativistic current sheets in self-consistent 
kinetic plasma simulations, identifying the generation of non-thermal particles 
and formation of a stable power-law shape in the particle energy distributions 
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f{'~f)d'~f oc ^^'"^drf. Depending on the dimension of the simulation, spectral in- 
dex from 2 (ID, attributed to a stochastic Fermi- type acceleration) to 3-4 
(recognized as a rather universal index of relativistic magnetic reconnection 
in previous 2D and 3D kinetic simulations, see Jaroschek at al. 2004, Zenitani 
and Hoshino 2005) were found. Lyubarsky and Liverts (2008) also studied 
the compression driven magnetic reconnection in the relativistic pair plasma, 
using 2.5D (i.e., 2D spatial, 3D velocity) simulations, finding that the spec- 
trum of particles was non-thermal, and a power law was produced. It seems 
a power-law distribution for the leptons inside the PWZ is then a plausible 
assumption. 

All in all, to find an a-priori dissipation solution for pulsars, and in particular, 
for the assumed pulsar in the LS 5039 which is the one we focus, is beyond the 
scope of this work (and actually, for the latter particular case, such solution is 
beyond what is by definition possible for a pulsar that we do not know exists) . 
In general, we note that an additional difficulty resides in the fact that the 
PWZ of pulsars in binaries may be subject to conditions others than those 
found in isolated pulsars. It is not implausible that close systems may trigger 
different phenomenology within the PWZ, ultimately affecting particle accel- 
eration there. Nevertheless, it is relevant for this paper to assume a particular 
particle's energy spectra with which we compute high energy processes in the 
PWZ, e.g. the up-Comptonization of the stellar field. We will assume two 
cases, a mono-energetic spectra -as a benchmark- (e.g., see Bogovalov & Aha- 
ronian 2000) and a generic power-law spectra (that could itself be subject to 
orbital variability). Both act as a phenomenological assumption in this paper, 
which goodness is to be assessed a posteriori, by comparison with data. 

As discussed, initial injection could come directly from the pulsar (the inter- 
acting particle population can of course be later affected by the equihbrium 
between this injected distribution and the losses to which it is subject, just as 
in the case of shock- provided electron primaries). The more compact the bi- 
nary is, the more these two settings (shock and pulsar injection, equilibrated by 
losses) are similar to each other. Given the directionality of the high Lorentz 
factor inverse Compton process, photons directed towards the observer are 
generally those coming from electrons moving in the same direction. Opaci- 
ties to processes such as inverse Compton and 77 absorption are high in close 
7-ray binaries, cascades can develop, and high-energy processes can already 
happen, as we explicitly show below, in the pre-shock region. 

In the model where pairs are injected as monoenergetic particles with the 
energy corresponding to the bulk Lorentz factor of the pulsar wind, they are 
frozen in the B-field. Under this assumption we neglect here the synchrotron 
losses, since there are none. 

When the injected pairs distribution is given by a power law spectrum the 
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r/R, r/R^ 



Fig. 2. Conditions upon the local magnetic field in the pulsar wind zone (PWZ). 
Left Panel: The black lines represent the maximum magnetic field, Bmax, for which 
IC losses dominate over synchrotron ones as a function of distance r/Rs from the 
massive star for 1 TeV electrons (thick curve) and 10 TeV (thin curve). The magnetic 
field in the PWZ, also as a function of radius from the massive star, is shown with 
a red line (for the pulsar at the apastron position) and with the blue line (for the 
pulsar at periastron). The system separation radius at periastron and apastron are 
also marked (gray lines). The magnetic field decays in the PWZ from its value 
close to the light cylinder. Right Panel: A different way of showing these conditions. 
PWZ magnetic field for different magnetization parameters {a = (labeled 
Bmax) 10-^ 10"^ (labeled B 

min) s-s a way of artificially introducing uncertainty 
in the real value, the magnetization parameter is expected to change within the 
PWZ (e.g., as in Crab). Maximum magnetic field for IC domination at periastron 
and apastron are again shown for different electron energies. The size of the star is 
represented by a shadow rectangle. The linear scale in the x-axis is in units of the 
light cylinder. The shock radius at periastron, and the separation at periastron and 
apastron are all shown. 

situation is more complex, since not all of them may be frozen in the PWZ 
field. The synchrotron cooling time is given by equation: 

tsyn = ^OOBa^ETlvS, (2) 



where Bq is the local magnetic field and is given in Gauss, and Et^v is electron 
energy given in TeV. For IC scattering, the timescale is instead given by 

U, = 7xl0^uJo^E^lyS (3) 



which gives good approximation for electron energy loss time for -E" ^ 1 TeV 
(e.g., Khangulyan et al. 2008). Timescales approach when 

E™'^^ ^ 0.24w°-^E^,V^ G. (4) 
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For a star with effective temperature Tg = 3.9 x 10^ K, the thermal field 
density at certain point at distance r from the massive star center is given by: 

cuo = AaT^/c X {Rs/2rf ^ 1.75 x 10^(i?,/2r)^ erg cm'l (5) 



Thus, the local magnetic field in the pulsar wind region is given by: 

5™"^ ^ 31.75{Rs/2r)E^^f G. (6) 



Applying this condition at periastron (r = 2.25Rs), the local magnetic field 
at the assumed injection place results in B^ax-per ~ 7 G. At apastron (r = 
4:.72Rs), it results in an stronger condition Bmax-apa ~ 3.4 G. That is, the 
magnetic field should be less than these values in order for IC to dominate 
over synchrotron losses at that particular position (the light cylinder) in the 
PWZ. Figure [2] shows this in detail. 



The magnetic field at the light cylinder distance Rlc is given by the dipole 
formula Blc = BoiRpsr / Rlc)^ ■ In the PWZ, the magnetic field is decreasing 
with distance as B(r) ~ ^a/{l + a) Blc{.Rlc l^)- For a millisecond pulsar 
we get Rlc ~ 5 x lO'' cm and Blc ~ 10^ G (assuming Bq = 10^^ G and 
P = 10 ms) up to Rlc ~ 5 x 10* cm and Blc ~ 8 x 10"^ G (assuming 
Bq = 10^^ G and P ~ 100 ms), where Rpsr ~ 10 km. The results for these 
different parameters are similar as there were obtained with fixed pulsar power 
in the model Lgd = 10^^ erg s~^ and they are related by the standard formula 
Lsd = B^R^g^c/ AR\(j- We also assume here that the magnetization parameter 
is 0" = 0.001, but have explored other values of this and other parameters as 
well, with similar results (see Figure [2]) . 



Given our results (see Fig. [2]) where we show the local magnetic field for which 
IC dominates and the magnetic field in the PWZ as a function of the distance 
from the light cylinder we can conclude that the injection for the model have to 
(generically) occur at some distance from the light cylinder, or/and, if closer to 
it, the synchrotron losses can be important. However, as the separation of the 
binary is ~ 10^^ cm, several orders of magnitude larger than the light cylinder 
(e.g., Rlc ~ 10* cm), the change of the injection place within e.g. ~ 1% — 5% 
already gives the initial injection distance at ~ (100 — 500)i?Lc- Thus, no 
significant effect in the PWZ photon spectra and lightcurve is produced. 
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2.3 Normalization of the relativistic particle power 



The fraction of the pulsar spin-down power ending in the e^e interacting 
pairs can then be written as: 



Assuming that the distance to the source is c? = 2.5 kpc, the normahzation 
factor for electrons traveling towards Earth is A = N^+e' / ^T^d? ■ The specific 
normalization factors in the expression of the injection rate N^+e- (E) will be 
given together with the results for two models in the corresponding sections 
below. In the models presented here, only a small fraction (~1%) of the pul- 
sar's Lsd ends up in relativistic leptons. This is consistent with ions carrying 
much of the wind energetics. In the case of mono-energetic lepton distribu- 
tion, where the energy of the primaries is fixed at Eq = 10 TeV, we have 
Ne+e-{E) oc S{E — Eq). In the case of a power-law in energy, that may be 
constant or vary along the orbit, we have N^+e-iE) oc E~"'\ 

2.4 On parameter interdependencies 

The rj expression represents the realization of a specific scenario fixed by differ- 
ent values of the pulsar and massive star parameters. For instance, assuming 
= 10^'^ ergs~^ and MgVs as given in the Table H] below, we get t] ^ 0.5 
for the periastron phase. But we would achieve the same value of rj with a 
smaller Lsd, say 10^^ ergs~^, if within uncertainties we adopt a smaller value 
of MgVs. Clearly, if only one of the quantities, Lsd or MgVs changes, the shock 
position moves, what results in a correspondingly smaller or larger PWZ in 
which cascading processes are set. A key parameter is then given by r^, instead 
of rj, since this is what defines the real size of the PWZ. The mild dependence 
of Ts on 7] makes possible to accommodate a large variation in the latter main- 
taining similar results. Fig. [3] shows this dependence for a relevant range of 



In addition, it is obvious that Lgd and (3 are linearly (inversely) related. We 
note that a fixed value of rj does not imply a specific value for Lgd as it is 
combined with parameters of the massive star wind, subject to uncertainty in 
their measurements for the specific case of binary treated (e.g., for LS 5039), 
if at all known. In the end, the same results can be obtained for different sets 
of parameters. There is a degeneracy between the shock position, which deter- 
mines the extent of the PWZ, and the injected power in relativistic electrons. 
For a smaller PWZ (e.g., if we assume a smaller Lsd for the same product of 




(7) 



13 



a. 



2 5 








2,0 








1,5 






apastron \/ 


1,0 








0,5 






periastron 


0,0 









10 10" 10 

L,j [ erg s"' ] 

Fig. 3. Dependence of the distance to the shock form the pulsar side, r^, as a function 
of the spin-down power, L^^^, for fixed values of star mass loss rate, M = lO"'^ Mq 
yr~^, and terminal star wind velocity, Voo = 2.4 x 10'^ kms~^. 

Ms^s) a larger amount of injected power compensates the reduced interact- 
ing region. We found that to get the similar results when the i] parameter is 
smaller, the /3 parameter have to be increased roughly by the same factor. 

2.5 Wind termination 

Very interesting in the context of LS 5039 model properties is the dependence 
of the distance from the pulsar to the termination shock in the direction to 
the observer as a function of the orbital phase. This is shown in Fig. |H We 
find that for both inclinations, the wind is unterminated for a specific range 
of phases along the orbit, i.e., the electron propagating in the direction of 
the observer would find no shock. The PWZ would always be limited in the 
observer's direction only if the inclination of the system is close to zero, i.e., 
the smaller the binary inclination the narrower the region of the wind non- 
termination viewed by the observer. In the discussed scenarios for LS 5039, 
the termination of the pulsar wind is limited to the phases ~ 0.36 — 0.92 
for inclination i = 60° and to ~ 0.45 — 0.89 for i = 30^, what coincides 
with the phases where the emission maximum is detected in the very high- 
energy photon range (see below the observed lightcurve obtained by H.E.S.S.). 
The unterminated wind is viewed by the observer at the phase range between 
apastron and INFC, while a strongly limited wind appears from periastron and 
SUPC. Note that the important differences in the range of phases in which the 
unterminated wind appears for distinct inclinations (e.g., the observer begins 
to see the unterminated pulsar wind at ~ 0.36 for i = 60° compared to ~ 0.45 
for i = 30^) produce a distinguishable feature of any model with a fixed orbital 
inclination angle. 

From the geometrical properties of the shock surface, there can also be spe- 



14 




Fig. 4. The distance from the pulsar to the termination shock in the direction to 
the observer (in units of stellar radius Rs), for the two different inclination angles, i, 
analyzed in this paper. INFC, SUPC, periastron, and apastron phases are marked. 
Additionally, the gray line shows the separation of the binary (also in units of Rg) 
as a function of phase along the orbit. 

cific phases for which the termination shock is directed edge-on to the ob- 
server. These are phases close to the non-terminated wind viewing conditions: 
slightly before and after those specific phases for which the wind becomes 
non-terminated, e.g., for the case of LS 5039 and i = 30°, (p ~ 0.4 and 0.93; 
whereas for i = 60°, it appears at ~ 0.32 and 0.93. Thus, we find that 
the condition for lepton propagation change significantly in a relatively short 
phase period, when the termination shock is getting further and closer to the 
pulsar (phases periods ~ 0.2 — 0.45 and ~ 0.9 — 0.96). In Table [1] we show 
the differences in these geometrical parameters for a few characteristic orbital 
phases. Note that even when they are essential to understand the formation 
of the lightcurve and phase dependent high energy photon spectra, all the 
above features are based on the simple approximation of the geometry of the 
termination shock. In a more realistic scenario, the transition between the ter- 
minated and unterminated wind, and its connection with orbital phases and 
inclination are expected to be more complicated yet. 



2.6 Opacities along the LS 5039 orbit 



The conditions for leptonic processes for this specific binary can be discussed 
based on optical depths to ICS and 77 absorption. The target photons for 
IC scattering of injected electrons and for e"'"e~ pair production for secondary 
photons are low energy photons of black-body spectrum with temperature 
Ts = 3.9 X 10^, which is the surface temperature of the massive star. This is 
an anisotropic field as the source of thermal photons differ from the place of 
injection of relativistic electrons, which for simplicity is assumed to be at the 
pulsar location for any given orbital phase. For fixed geometry parameters, 
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Table 1 

Geometrical properties for specific phases afong the orbit, pi and p2 stand for 
the phases at which the angle to the observer is 90°; ri and r2 stand for the phases 
between which the shock in the observer direction is non-terminated. The separation 
d and the shock distance are given in units of star radii Rg. 
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Fig. 5. Opacities to ICS for electrons and pair production for photons as a function 
of electron/photon energy. Opacities were calculated up to infinity for injection at 
the periastron (left panel) and apastron (right panel) pulsar distance for different 
angles of propagation with respect to the massive star. 

the optical depths change with the separation of the binary (in general with 
the distance to the massive star), the angle of injection (the direction of prop- 
agation with respect to the massive star), and the energy of injected particle 
(the electrons or photons for 7 absorption). 

To have a first handle on opacities, we have calculated the optical depths 
adopting the binary separation at periastron and apastron. In case of the 
optical depths for photons we had to do an additional simplifying assumption 
to allow for a direct comparison with the optical depth for electrons. In the 
discussed model of this paper, photons are secondary particles and they do 
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Fig. 6. Opacities to ICS for electrons and pair production for photons as a function 
of the injection angle, a, in the LS 5039 system, for different initial energies of the 
particles (left panel E = 1 TeV, right panel E = 10 TeV) and distance to the massive 
star (including periastron and apastron separation). Opacities were calculated up to 
infinity. The injection angle is provided in the plot and for pair production process 
they are given in the same order as for ICS. 

not have the same place of injection as electrons, but because the path of 
electrons up to their interaction is much smaller than the system separation, 
we assume now the same place of injection for photons and electrons: for a 
first generation of photons, this is sufficient for comparison. 

Results are shown in Fig. [5] and the needed formulae for its computation are 
given in the Appendix. The optical depths for ICS can be as high as ~ few 
100 for significant fractions of the orbit, and are still above unity for energies 
close to 10 TeV. With respect to the injection angle, the opacities are above 
~ 1 for all propagation angles, except the outward directions at the apastron 
separation. (Recall that outward -inward- directions refer to a-angles close to 
0° and 180°, respectively, as shown in Fig.[l]). In the case of 77 pair production, 
the interactions are limited to specific energies of the photons and are strongly 
angle-dependent. The most favorable case for pair creation is represented by 
those photons with energy in the range 0.1 to few TeV, propagating at > 50°. 
The highest optical depths to pair production are found for directions tangent 
to the massive star surface: for propagation directly towards the massive star 
pair production processes are limited by the star surface itself (for periastron 
the tangent angle is at ~ 154°; for apastron, it is given at an angle of ~ 168°). 
This is better shown in Fig. iGlHI 

The opacities for pair production are characterized by maxima which change 
with respect to the angle of photon injection. For inward directions (in general 



The range of parameters depicted in this figure is just to emphasize that there 
could also be propagation and interactions outside the innermost parts of the system, 
in general, although of course with a reduced opacity. The actual range of angles 
towards the observer allowed in a particular system depend on the orbit inclination. 
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Fig. 7. Opacities to ICS for electrons and pair production for photons as a function 
of the orbital phase, in the LS 5039 system, for different inclination of the binary 
orbit (left panel i = 30'', right panel i = 60'^) and energy of injected particle. 
Opacities are calculated for specific phases up to the termination shock in direction 
to the observer. 

for angles larger then 90°), the peak energy is higher, ~ 10^ GeV— 3 x 10^ GeV, 
than for outward directions, where we find ~ 3 x 10^ GeV— 5 TeV. For energies 
up to the peak for pair production, the IC process dominate for the same angle 
of particle injection, but for higher energies the opacities for 7-absorption 
become slightly higher than for ICS. 

Cascading processes are effective when the interaction path for electrons and 
photons are short enough so that a cascade can be initiated already inside the 
pulsar wind zone. As we can see in Fig. [5], the opacities at high-energy range 
(Klein-Nishina) are comparable to the maximal opacities for 77 absorption. 
This means similar probability of interaction for both photons and electrons. If 
the electrons scatter in the Klein-Nishina regime, the photon produced, with 
comparable energy to that carried by the initiating electron, is likely to be 
absorbed, and a e+e~ cascade is produced. 

The dependence of the optical depths upon the orbital phase for specific pa- 
rameters of the LS 5039 is shown in Fig. [3 These opacities are calculated up to 
the termination shock in the observer direction. The presence of the shock (at 
a distance r^) limits the optical depths. Since this parameter is highly variable 
along the orbit (see Fig. H]), its influence on the optical depth values is not 
minor. As the angle to the observer, which defines the primary injection angle, 
vary in the range {90 — i,90 + i) (for INFC and SUPC, respectively), we found 
that there is a range of orbital phases for which the PWZ is non-terminated. 
In that case the cascading process -which develops linearly- is followed up to 
electron's complete cooling (defined by the energy -Emm = 0.5 GeV). The op- 
tical depths in the non-terminated pulsar wind are thus comparable to those 
presented in Fig. [5l taking into account the differences between the separation 
and angle of injection. 



18 





Fig. 8. Comparison of the opacities to ICS for electrons and pair production for 
photons along the orbit of LS 5039, for fixed particle energy but different inclination 
angles (i = 30°, z = 60°). 



Apart from the condition for electron cooling, in the case of the terminated 
wind the cascade is followed up to the moment when the electron reaches 
the shock region, whatever happens first. As can be seen from comparison of 
Fig. [5] and [71 the opacities are significantly smaller when the propagation is 
limited. For instance, as a comparison we can choose the electron injection 
at periastron, the angle aobs = 110° (130° for i = 60°), and fix the energy 
to ~ 1 TeV. Then, the opacities along the orbit up to the terminated shock 
are ~ 2 — 3 while for the unterminated processes the opacity are ~ 10 — 20. 
A direct comparison for specific energies of injected particles and the two 
inclinations considered is shown in Fig. [HI Moreover, the opacities for both 
processes decreases along the propagation path, see Fig. [HI as the angle to the 
observer also decreases. This indicate that even if the photons propagate very 
close to the massive star, the probability for absorption can be much smaller 
than that at its injection place, and finally become less then 1 at the distance 
from the pulsar ~ 3 — 4 i?^, for an energy of 10 TeV and 1 TeV respectively. 
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Fig. 9. Opacities to ICS for electrons and pair production for photons as a function 
of the distance from the injection place along the propagation path from the pulsar 
side, dr, in the LS 5039 system. Opacities are calculated for two values of the 
particles initial energies, E = 1 TeV (left panel), and £^ = 10 TeV (right panel), 
and different angles of injection. 

Table 2 

Model parameters for the example of mono-energetic injection considered 



Meaning 


Symbol 


Adopted value 


Spin-down power of assumed pulsar 






Fraction of Lsd in relativistic electrons 


P 


10-2 


The energy of injected pairs 


Eq 


10 TeV 



3 Mono-energetic electrons in the PWZ 

Following the normalization formula given by Eq. ([7]), for the primary lep- 
ton energy E = TeV and the nominal value of pulsar spin-down power 
Lsii = 10^^ergs-\ and assuming a distance to the source d = 2.5 kpc we get 
ATg+g- = px 6.24 X 10^^ 5{E - lOTeV) s"^. The normalization factor for elec- 
trons traveling towards Earth is A = Ne/4:7rd'^ = /? x 8.83 x 10"^° s^^ cm^'^. As 
a first normalization of the simulation results below, we have fixed (3 = 10^^. 

3.1 Lightcurve and broad-phase spectra 

Based on the simulations for the specific phases along the orbit, the photon 
spectra and lightcurves in different energy ranges (10 — 100 GeV, 10 — 10^ GeV, 
and above 1 TeV) were calculated for two assumed inclinations of the orbit 
of the system, i = 30° and i = 60°. The highest energy range corresponds 
to the data presented by H.E.S.S. (Aharonian et. al. 2006). To allow for a 
comparison, the photon spectra were calculated also for a primary electron 
energy of 1 TeV for specific orbital phases (periastron, apastron, SUPC and 
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Fig. 10. VHE photon spectra for specific phases, at periastron, apastron (thin lines), 
INFC, and SUPC (bold lines), for two inclination angles, i = 30° (top) and i = 60° 
(bottom). The spectra were calculated for two different primary energy of electrons, 
ii^ = 10 TeV (left) and E = 1 TeV (right). The few free parameters involved in 
the mono-energetic model and their assumed values are given in Table [2j For direct 
comparison, the same normalization factor was used in case of each injection energy. 
Discussion on the non- uniqueness of model parameters is given in Section 2. 

INFC). 

The photon spectra (SED) are presented in Fig. [101 AH the spectra for primary 
electron energy E = 10 TeV are hard, presenting a photon index < 2, and 
they have significantly higher flux at the highest energy range, close to the 
initial primaries' ~ 10 TeV, as the ICS of primary pairs mainly occur in 
Klein-Nishina range. Below this peak the spectra are well represent by the 
power-law with photon index close to ~ —1, albeit one can also notice the 
change in the photon index from SUPC to INFC (harder spectrum) and the 
notably different energy cutoffs. These features are discussed below. 

While for phases around SUPC the PWZ is limited in the direction to the 
observer, for the opposite phases the PWZ is unterminated. On the other hand, 
for INFC the angle to the observer is at its minimum (and also the separation 
of the system is larger) what causes the decreasing of optical depth to both 
processes considered. In contrary to the INFC phases, those at SUPC present 
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Fig. 11. Theoretical lightcurves for mono-energetic injection and two different incli- 
nation angles. As before: i=60'', solid; i=30'', dot-dashed. 

strong absorption features at an energy range from ~ 0.1 to few TeV, which 
cause the spectra to be cut at lower energies. This is due to the dependence of 
the optical depth to 77 absorption on the energy of photon and the influence 
of geometry in defining the photon path towards the observer. The larger the 
angle to the observer, the stronger the effects of absorption are, what can 
be seen from comparison of the results for two inclination angles (SUPC and 
periastron phases). Notice, that for i = 60° the angle to the observer attains 
the largest value of the discussed examples, aobs = 150". The absorption for 
SUPC phases takes place already in the PWZ, what can initiate the cascading 
processes. Photons which pass through the shock region can be also absorbed 
in the MSWZ. The opacities for 77 absorption also strongly depend on photon 
energy (see Fig. [8]). High energy photons, with energy from few GeV up to 
few TeV, are the most likely to be absorbed, as the opacities decrease with 
energy for all shown angles of the photon injection. Indeed, the number of the 
photons absorbed in MSWZ is about ~ 10 — 20 % of the number of photons 
reaching the shock region and the effect is significant in the final spectra. 
On the other hand, the cascading in PWZ produce lower energy leptons and 
photons (in the cooling process of secondary pairs) and this causes a higher 
photon flux at energies up to few 10 GeV. The processes in PWZ are limited 
by the presence of the shock, what influence the photon production rate. For 
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INFC phases the absorption of 7-rays is minor, as the photons propagate 
outwards of the massive star; once produced, most photons can escape from 
the system (also because of the threshold to pair production). Higher energy 
photons are produced mostly in the first IC interaction of the primary pairs, 
while further electron cooling, not limited by the termination shock, supply 
the spectrum in lower energy photons. This flux is not as high as in the case 
of SUPC phases, where most of the cascading takes place. 



Similar dependencies, both for INFC and SUPC phases, are present in the 
spectra obtained from the simulations for primary energy of pairs E = 1 TeV. 
As the optical depths to IC scattering are higher in that case, the processes 
are more efficient and the number of produced photons are higher. To give an 
example, most photons of energy > 20-30 GeV at SUPC-periastron phase are 
absorbed in MSWZ. 



Comparing the photon spectra produced at SUPC and INFC, an anticorrela- 
tion between GeV and TeV photon fluxes is evident, at least comparing the 
spectra at energies below and above ~ 10 GeV. All these effects play a role 
in the formation of the 7-ray lightcurve, shown in Fig. [TTJ The lightcurve in 
the highest energy range (> 1 TeV) has a broad minimum around SUPC, 
0.96 < < 0.25. From the comparison of the position of the termination 
shock along the orbit, the opacities to ICS and 77 absorption we see that 
this minimum is formed despite being at phases with the highest opacities to 
both processes considered and so having an effective cascading in the PWZ 
region. This minimum is mainly due to the absorption of the photons which 
get through the shock and propagate into the massive star wind region. The 
broad maximum in the high-energy lightcurve is formed at opposite phases 
around INFC, 0.4 < (p < 0.9. The maximum for higher inclination is also 
characterized by a local minimum close to the INFC phase, what is the re- 
sult of the IC opacity dependence on the propagation angle. The range of 
INFC phases corresponds to the unterminated pulsar wind in the observer 
direction. The opacities for both inclinations have local peaks in this range, 
what was discussed in the previous paragraph. The local minimum in photon 
flux within the INFC range reflects a similar behavior of the opacities. So, 
as far as the propagation of the particles in the PWZ is unterminated, the 
high-energy lightcurve formation is in agreement with the dependence of the 
optical depths. Additionally, we can also see that the flrst local peak in this 
lightcurve is formed earlier in phase for inclination i = 60° (0 ~ 0.35), what is 
also in agreement with the dependence shown by the distance to the termina- 
tion shock. Note that the second peak in the broad TeV lightcurve maximum, 
at ~ 0.85, is higher than the flrst one, at ~ 0.4 what is the result of a change 
in the separation of the system together with the angle to the observer. 



23 



3.2 Peaks and dips in the lightcurve 



The geometrical conditions for leptonic processes change significantly along 
the orbit, and they are more efficient for phases around SUPC. However, 
these processes, as we can already see from the opacities dependence on the 
orbital phases, is limited by the presence of the termination shock. This causes 
that the best conditions for very high-energy photon production, and finally 
escaping from the system, occur for an specific combination of the angle to 
the observer, the separation of the system, and the distance to the shock from 
the pulsar side. From Figs. Hand H] we can see that this happens at the phases 
~ 0.3 — 0.4 and ~ 0.9 what refiects in the TeV lightcurve. At phases close to 
INFC photons are produced mainly in the primary e~^e~ pairs cooling when 
the propagating electron undergo frequent scatterings. Together with the fact 
that there is no efficient absorption of photons once they are produced, they 
finally escape from the system, what yields to the broad maximum in the TeV 
lightcurve. On the other hand, at phases around SUPC, high-energy photons 
are absorbed when propagating through the system in the MSWZ, what causes 
a dip in the lightcurve. For these phases many more lower energy photons are 
produced in cascades in the PWZ, what yields to a maximum in the GeV 
lightcurve for SUPC, anticorrelated with the behavior at TeV energies. 



4 A power-law electron distribution 

For further exploration of the 7-ray production model we set a new assump- 
tion: the energy distribution of the interacting e^e~ pairs is given by a power- 
law spectrum. We will additionally assume that the power-law may be con- 
stant or vary along the orbit. Motivated by the different observational behavior 
found, we have assumed that two different spectral indices correspond to the 
two broad orbital intervals proposed by HESS (Aharonian et. al., 2006). For di- 
rect comparison we have specified the interval around the inferior conjunction 
(with the apastron phase): 0.45 < (j) < 0.90, and around superior conjunction 
(including the periastron phase): < 45 and (p > 0.90 as being bathed by 
different electron distributions. The results for power-laws distribution were 
already summarily presented in our earlier work Sierpowska-Bartosik and Tor- 
res (2007a,b) and we refer to these works for further details. The agreement in 
both spectra and lightcurve is notable (particularly for the case of a variable 
lepton spectrum along the orbit), as can be seen in Figs. [T^ and [TBI discussed 
below. 

As we could have already noticed from Fig. [5l the optical depth for high 
energy electrons is below unity. In that case, part of the initial electrons will 
be interacting in the PWZ less efficiently and finally will reach the shock 
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Table 3 

Model parameters for interacting electrons described by power-laws 



Meaning 


Symbol 


Adopted value 


Spin-down power of assumed pulsar 




Xo37 gj.g s^^ 


VHE cutoff of the injection spectra 


E 

^max 


50 TeV 


Constant lepton spectrum along the orbit 


Fraction of Lgd in leptons 


P 


10-2 


Slope of the power-law 


Fe 


-2.0 


Variable lepton spectrum along the orbit 


Fraction of Lsd in leptons at INFC interval 


/? 


8.0 X 10-3 


Slope of the power-law at INFC interval 




-1.9 


Fraction of Lg^ in leptons at SUPC interval 




2.4 X 10-2 


Slope of the power-law at SUPC interval 




-2.4 




E [GeV] E [GeV] 



Fig. 12. Comparison of the spectra of electrons injected in the PWZ with the cor- 
responding spectrum of electrons which reach the shock region (SUPC, periastron) 
or leave the innermost part of the PWZ after propagation in it (INFC, apastron) 
in case of non-terminated pulsar wind. Left: the electron spectrum injected at the 
SUPC and periastron (PER) with the initial index Fg = —2.4 (dot-dashed line) and 
spectrum for electrons reaching the shock: at SUPC (solid line) and PER (dashed 
line). Right: the electron spectrum injected at the INFC and apastron (APA) with 
the initial index Fg = —1.9 and spectrum for electrons after interaction in the local 
radiation field (the shock for this phases is not terminated in the direction to the 
observer): at INFC (solid line) and APA (dashed). 

region (we remind that in the direction of the observer, there is not always a 
shock in the electron's propagation). In Fig. [12] we show the spectra of initial 
electron distribution and corresponding electron spectra after propagation in 
the PWZ. The electrons which reach the termination shock are isotropised 
there. The termination shock shape is specific for each phase and it is limited 
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Fig. 13. The spectra of photons escaping from the system produced by spectrum of 
electrons injected at separation distance (sohd hues) and inside the PWZ (dashed 
hues) at SUPC (left figure) and INFC (right). Left: the photon spectrum produced 
at SUPC {(j) ~ 0.06) by electrons injected inside PWZ, at Vinit = (the spectrum 

index Tg = —2.4). Right: the photon spectrum produced at INFC {(p ~ 0.72) by 
electrons injected inside PWZ, at Vinit — 

D- (Te = -1.9) as the shock at INFC 
is not terminated in the direction to the observer. 

in space. The electrons at the shock, locally re-accelerated, become the initial 
spectra for the next generation of photons (not only radiative but adiabatic 
losses have to be taken into account). 

The MSWZ can play a role (although we believe it will not be too impor- 
tant in such close binaries like LS 5039, because of the high opacities already 
encountered in the PWZ along most of the orbit and the loss of direction- 
ality -consequently, of random photon emission- of the electron population). 
An estimation of the contribution of the cascades in MSW is not trivial, espe- 
cially from the normalization point of view, as the magnetic filed in this region 
causes isotropisation of produced photons. In addition, the magnetic field of 
the massive star is ordered, what could cause additional effects e.g. focusing 
the propagating electrons in some regions close to the massive star (where 
the magnetic field is dipolar) as was shown in Sierpowska & Bednarek (2005). 
The impact of the MSWZ could make the need for a change in the injection 
index of relativistic particles less severe. A 3D cascading code is needed for 
such estimates, and we expect to report on that in the future. 

The model assumes that the initial electron are injected in the vicinity of the 
pulsar light cylinder. This is actually an assumption which we have investi- 
gated further. Indeed, we have investigate also the changes in the produced 
photon spectra if injection take place at a further distance inside the PWZ. 
Results are given in Fig. [131 For phases where the wind is terminated in 
the direction to the observer (e.g., SUPC) the new injected radius was fixed 
to Tinit = D — rsh/2 from the massive star, where D is separation at given 
phase (i.e., at the middle of the PWZ). For phases such as INFC, where the 
shock is unterminated, the injection place was shifted by the radius of the star 
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Fig. 14. Very high-energy spectra of LS 5039 around INFC and SUPC, together 
with the theoretical predictions in equal phase intervals for power-laws electron 
distribution and two different inclination angles. The free parameters involved in 
our model and their assumed values are given in Table El After Sierpowska-Bartosik 
& Torres (2008). 
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Fig. 15. H.E.S.S. run-by-run folded (ephemeris of Casares et al. 2005) observations 
(each point corresponds to 28 min of observations) of LS 5039 with the results of 
the theoretical model for power-law distribution and two different inclination angles 
for 1=30*^ (left) and i=60^ (right). Light (green) lines stand for results obtained with 
a constant interacting lepton spectrum along the orbit, whereas dark lines (black) 
correspond to a variable spectrum. After Sierpowska-Bartosik &: Torres (2008). 

Tinit = D — Rg, that also corresponds to the half-distance of the separation 
between the pulsar and shock r^/j for this phase, if movement is in the direc- 
tion to the massive star. We can see that the produced photon spectra do not 
differ significantly in this two scenarios, making our results stable. 



4-1 Lightcurve and broad-phase spectra 



To construct the lightcurve, spectra for over 20 orbital phases were calculated 
to cover the whole orbit of LS 5039. The specific averaged spectra were ob- 
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tained based on the same orbital intervals as presented in H.E.S.S. data. The 
averaged spectra were obtained summing up individual contributions from or- 
bital spectra in given phases, each with a weight {5t/T) corresponding to the 
fraction of orbital time that the system spends in the corresponding phase 
bin. Then, comparing the observational and theoretical averaged spectra, the 
parameter (3 was estimated. Both for constant and variable injection model, 
we get that the fraction of the spin-down power in the primary leptons has 
to be at the level of ~ 1 %. With this parameter in hand each of the single 
spectra can be equally normalized. 

It is worth noticing how well these lightcurves compare with those in the work 
by Bednarek (2007), at least for some of the specific phases considered by him. 
Bednarek also included cascading in his simulations, and the physical input 
of his model (although in the case of a microquasar scenario) is similar to 
ours. As a result, the anti-correlation phenomena (from GeV to TeV energies) 
is also a result of his work. The spectrum along the orbit with respect to 
H.E.S.S. datapoints and the possible short-timescale variability (see below) 
was not provided by Bednarek, so that a comparison with these results is not 
possible. A detailed differentiation between these two models is a key input 
for distinguishing (microquasar or 7-ray binary) scenarios, even when some 
assumptions are intrinsic to each of the models, isolating the contribution of 
an equal physical input can help decide on what object constitute the system 
LS 5039. As an example: Bednarek did find in his model that a fixed inclination 
angle {i = 60°) was needed in order to reproduce the shape of the H.E.S.S. 
hghtcurve results, whereas in our case, as we see in Figure [13 the influence of 
inclination is minor. 

For completeness, we mention that as noted above, H.E.S.S. has also provided 
the evolution of the normalization and slope of a power-law fit to the 0.2-5 TeV 
data in 0.1 phase-binning along the orbit (Aharonian et al. 2006). The use of a 
power-law fit was limited by low statistics in such shorter sub-orbital intervals, 
i.e., higher-order functional fittings such as a power-law with exponential cutoff 
were reported to provide a no better fit and were not justified. To directly 
compare with these results, we have applied the same approach to treat the 
model predictions, i.e., we fit a power-law in the same energy range and phase 
binning. We show here this comparison in Fig. [12] (taken from Sierpowska- 
Bartosik and Torres 2008) in the case of a variable lepton distribution along 
the orbit. We find a rather good agreement between model predictions and 
data. Results for constant lepton distribution along the orbit can be seen in 
Sierpowska-Bartosik and Torres (2007a). 

In Fig. [TUthe SED in H.E.S.S. energy range are shown for both, constant and 
variable lepton distributions along the orbit, and two inclination of the system 
i = 30° and i = 60°. It was shown in the previous Section that in case of the 
mono-energetic injection, there are photon flux differences in INFC phase due 
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Fig. 16. Shaded areas in the left (right) panel show the change in the normalization 
(photon index) of a power-law photon spectra fitted to the theoretical prediction 
for each of the 0.1 bins of phase in the case of a variable lepton distribution along 
the orbit. The two different colors of the shading stand for the two inclination 
angles considered. The size of the shading gives account of the error in the fitting 
parameters. Data points represent the H.E.S.S. results for a similar procedure: a 
power-law fit to the observational spectra obtained in the same phase binning. 
From Sierpowska-Bartosik & Torres (2008). 

to the change in the angle to the observer. In the power-law distribution model 
this effect is less significant due to the contributions of different energies, and 
the dependence on energy of the processes involved. Based for instance on the 
injection of constant spectrum with slope Tg = —2.0 we can see that there is 
no significant differences involving the inclination. 

Exploring the more evolved model of the variable injection we see that the 
steeper the primary spectrum is, the higher the photon flux at lower energy 
range (SED for lower energies were shown by Sierpowska-Bartosik and Torres 
2008). For constant lepton distribution along the orbit, the differences in pho- 
ton flux below 100 GeV are not so significant, but the difference gets important 
in the case of a variable lepton distribution (the difference between both mod- 
els in the lower energy range is about one order of magnitude). Overall, the 
variable lepton distribution provides a better agreement with all data (note 
that it has only two extra free parameters when compared with the constant 
distribution case, see Table [3l but matches more than 10 data points that were 
missed in the previous case). In particular, it is worth noticing that there is 
good agreement with the H.E.S.S. spectra at an energy ~ 200 GeV where both 
SUPC and INFC spectra coincide and above which the photon flux for INFC 
is higher then the flux for SUPC dominating in the lower energy range; i.e., the 
energy where the anticorrelation begins (see Sierpowska-Bartosik and Torres 
2008). It is also interesting to remind that absorption alone would produce 
strict modulation (zero flux) in the energy range 0.2 to 2 TeV whereas the ob- 
servations show that the flux at ~ 0.2 TeV is stable. Additional processes, i.e. 
cascading, must be considered to explain the spectral modulation. Our model 
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Fig. 17. Evolution of individual spectrum in the case of variable lepton distribution 
along the orbit, at individual phase bins, from 0.1 to 0.9. The shadow represents 
the GLAST energy coverage; the rest of the span of the x-axis can be observed by 
ground-based facilities. 

have these processes consistently included and the lightcurve details arise then 
as the interplay of the absorption of 7-rays with the cascading process, in the 
framework of a varying geometry along the orbit of the system. 



4-2 Testing with future data 



Apart from possible testing with GLAST (at the level of lightcurve, spectra, 
hardness ratios, and differentiation between constant and variable electron 
distribution, see Sierpowska-Bartosik and Torres 2008), we can provide further 
possible tests at high 7-ray energies. It was already said that power-laws do not 
always present the best fit to the specific spectra along the orbit. Even when 
fitting such power-laws to the theoretical predictions provides agreement with 
data (see Fig. dn]), observations with larger statistics (with H.E.S.S., H.E.S.S. 
II, or CTA) could directly test the model in specific phases. A model failure 
in specific phases would allow further illumination about the physics of the 
system. Figure [T7I shows the evolution of individual spectrum in the best fitting 
case of variable lepton distribution along the orbit, at individual phase bins. 
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from 0.1 to 0.9, for testing with future quality of data. 



5 Concluding remeirks 

We presented the details of a theoretical model for the high-energy emission 
from close 7-ray binaries, and applied it to the particular case of LS 5039. 
The model assumes a pulsar scenario, where either the pulsar or a close-to-the- 
pulsar shock injects Icptons that after being reprocessed by losses to constitute 
a steady population, are assumed to interact with the target photon field 
provided by the companion star within the PWZ. The model accounts for the 
highly variable system geometry with respect to the observer, and radiative 
processes; essentially, anisotropic Klein-Nishina ICS and 77 absorption, put 
together with a Monte Carlo computation of cascading. The formation of 
lightcurve and spectra in this model was discussed in detail for the case where 
the interacting leptons are assumed mono-energetic and described by power- 
laws. 

Comparing the interacting models of mono-energetic leptons and power-law 
distributions we can see similar dependencies for the spectral changes along 
the orbit and the GeV to TeV lightcurves. Thus, the case of a mono-energetic 
population is useful to understand some of the aspects regarding the forma- 
tion of the observational features, although it does not match observational 
data. For the mono-energetic case, we have shown the spectra produced at 
characteristic phases (INFC, SUPC, periastron, apastron) for primary ener- 
gies of 1 TeV and 10 TeV. In the power-law distribution model, the specific 
features of the photon spectra and the lightcurves produced for an specific 
primary electron energy overlap. Wc can still notice the absorption features 
in the spectra produced at SUPC. We have discussed effects solely based on 
the optical depths and the general geometrical dependence along the orbit, as 
well as on the presence of the shock which terminates the pulsar wind in the 
direction to the observer at SUPC phases. 

A power-law lepton distribution interacting in the PWZ describes very well 
the phenomenology found in the LS 5039 system at all timescales, both flux 
and spectrum- wise, even at the shortest timescales measured. This latter re- 
sult is unexpected: wc find that there is nothing a priori in the model that 
allows one to predict that when broad phase spectra data (INFC and SUPC) 
are reproduced so will be the data at the individual and much shorter phase- 
binning, less with such a good agreement. This result point perhaps to some 
reliability of the model, at least in its essential ingredients: geometry, cascad- 
ing, interacting electron population. 

However, this model certainly has room for improvement. We emphasize here 
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that we do not have an a priori model for the interacting lepton population 
itself, although we have discussed the research on dissipation processes in the 
PWZ which may give raise to such distributions if it results from pulsar injec- 
tion. In any case, the assumption of power-laws is an approximation to a more 
complex scenario where the real interacting lepton population is the result of 
a full escape-loss equation. In addition, we are not considering yet the multi- 
wavelength emission at lower energies, since wc left out of our description the 
synchrotron emission of electrons accelerated at the shock and the morphology 
of the shock along the orbit, what we expect to discuss elsewhere. Other than 
system scalings that are fixed by multiwavelength observations, the model is 
based on just a handful of free parameters, and it is subject to tests at high 
and very high-energy 7-ray observations with both GLAST (described in more 
detail in Sierpowska-Bartosik & Torres 2008) and future samples of data at 
higher energies, where more statistics at finer phase bins can determine better 
the spectral evolution along the orbit of this interesting system. 

We acknowledge extended use of lEEC-CSIC parallel computers cluster. We 
acknowledge W. Bednarek for discussions, and the Referee for useful com- 
ments. This work was supported by grants AYA 2006-00530 and CSIC-PIE 
2007501029. 
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Appendix: Numerical implementation and formulae 

This Appendix introduces further essential details concerning geometry and 
formulae for the implemented process of Inverse Compton scattering and e^e~ 
pair production in the anisotropic radiation field of the massive star. 

Monte- Carlo implementation for the cascading process 

A Monte-Carlo procedure is applied to calculate the place of electron inter- 
action, Xe, and the energy of the resulting photon in the ICS, E^; as well 
as the place of photon interaction, Xp, and the energy of the produced elec- 
tron/positron, (see, e.g., Bednarek 1997). The following summarizes the 
procedure. 

When the primary electron is injected, the computational procedure for IC 
scattering is invoked first in order to get the initial place of interaction (the 
radial distance from the injection position) and the energy of the up-scattered 
photon if the interaction takes place. For this same electron, the procedure is 
repeated afterwards up to the moment of electron cooling or when reaching 
the shock region. The photons produced along the electron path switch on a 
parallel computational procedure for pair production. This in turn gives the 
place of photon absorption and the energy of the e~^e~ pair if the process 
occurs. Because the photons can get trough the termination shock if not ab- 
sorbed in the PWZ, the numerical procedure is not limited to the termination 
shock radius. If a e'^e" pair is created within the PWZ, then the procedure 
for IC scattering is initiated for it, from which a next generation of photons 
can be produced. When photons cross the shock, information about them is 
separately saved in order to get account of the level of fiux absorption in the 
MSWZ. The emission spectra are produced from photons which finally escape 
from the binary, so the photons absorbed in the MSWZ are not included in it. 

The place of lepton interaction in IC scattering, that is, the production of 
a photon at Xei, after being injected at a distance Xi from the star, at an 
angle a (we discuss further details of geometry in the next Section and in the 
Appendix, especially see Fig. [23D, is calculated from an inverse method. For 
an specific random number Pi in the range (0, 1), the interaction place Xei is 
given by the formula: 



where XjQg{Ee,Xi,a,Xe) is the rate of electron interaction to IC process (see 
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Fig. 18. Left: The distribution of the place of electron interaction due to IC scaterring 
in the anisotropic radiation field of a massive star. The initial electron energy is 
here assumed as E = 1 TeV, and the place of its injection is given by the distance 
to the massive star in the LS 5039 system, ds = 2.25 Rs, being injected at the 
angle a = 150° with respect to the massive star. The number of injected electrons 
(interactions) is N = 3360. Right: Distribution for the scatterred photon energies 
for the same simulations. 



Eq. [95]of the Appendix) . From this it follows: — In Pi 



ICS 



where the integration is over the propagation path of the electron, Xe- Knowing 
vg, Xi, a, Xe), we keep integrating forward until we find numerical equal- 



ity of the integral result with the random-generated quantity — In Pi, what 
defines the position of interaction x^i- For a large number of simulations, the 
distribution of interaction distances is presented in Fig. f|T8| left panel). 



Once we get the place of electron interaction, Xgi, simulated in the way de- 
scribed, we get the energy of the scattered photon from an inverse cumulative 
distribution. The energy of the photon E'^ produced in IC is then given by the 
spectrum of photons after scattering dN /dE^dt (see Appendix, Eq. |H5]). For 
a random number P2 we define: 



P, = 



dE^ 

dE^dt ^ 



dE^dt ^ 



(9) 



where E^"-^ is the maximal scattered photon energy (see Appendix Eq. [72]) . 
while the energy E^ is the Monte-Carlo result for the simulated photon energy. 
Note that the denominator of the latter expression is the normalization needed 
for the Monte-Carlo association to succeed. The statistics for the energies of 
up-scattered photons is shown in Fig. [18], right panel. All relevant formulae 
and more details of implementation are given in the Appendix. 

In a similar way we randomize by Monte-Carlo the needed magnitudes for 77 
absorption. The probability of interaction for a photon with energy E^ at a 
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given distance, say Xps, from an injection place (i.e., provided Xi and a are 
known) during the propagation in an anisotropic radiation field is given by 
expression: 



exp 



Xp3 

- / A~i(E^,Xj, 



77 



(10) 



From this it follows, /q '''' X~^{Ey, Xi, a, x^) dx^ = — In P3, thus, similarly to the 
process just described, for the random number P3 we can get the specific place 
of interaction Xps. The results of such Monte-Carlo simulations are presented 
in Fig. [ini left panel. 

The energy of the lepton (e^e^) produced in the photon absorption process 
(with the photon having an energy E^) is calculated again from an inverse 
function. The latter is obtained from the integration of the e^e~ spectra pro- 
duced in the process: (see Appendix, Eq. |67|) : 



Pa = 



cEl dW{E~f,xo,a,Xp) ^ 

r^r"" dW{E-,,xo,a,Xp) ^ jf 
Jo.bE^ dEedx-y "-^e 



where P4 is a new random number. The integration is over the electron en- 
ergy Eg. For normalizing, as the produced lepton spectra are symmetric with 
respect to the energy E^ = E^/2, the lower limit of this integral is fixed to 
this latter energy, while the upper limit is equal to maximum energy available 
in this process E'^"'^ = E^ — me(?- From the energy of the electron, E^, the 
energy of the associated positron is also obtained as E^' = E^ — E^. Again, 
note that the denominator of the latter expression is the normalization needed 
for the Monte-Carlo association to succeed. 



Checks for the Monte- Carlo distributions 



Here we discuss the rightness of the simulated random distribution of electrons 
and photons resulting from the cascading process. In Fig. [20]we show the com- 
parison between the event statistics, i.e., N{x < Xe)/Ntotai with Ntotai being 
the total number of simulations run in these examples, and the analytical com- 
puted probability of interaction 1 — exp(— a;/A), where A is, correspondingly, 
the one corresponding to 77 absorption and ICS. We see total agreement of 
the Monte-Carlo and analytical probabilities, i.e., whereas the position of in- 
teraction of a single photon or electron, from which the subsequent cascading 
process is followed, is obtained through Monte-Carlo and thus is random, the 
overall distribution maintains the shape provided by the physical scenario: 
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Fig. 19. Left: Distribution of the position of the 7-photon interaction due to absorp- 
tion in the anisotropic radiation field. The initial photon energy in this example is 
assumed as E = 1 TeV and the place of its injection is given by the distance to 
the massive star in the LS 5039 system, ds = 2.25 Rs at the angle a = 150° with 
respect to the massive star. The number of injected photons (interaction) is N = 
3358. Right: Distribution of the produced electron energy for the same simulations. 
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Fig. 20. Comparison between the event statistics, i.e., N(x < Xe)/Ntotai with Nt^tai 
being the total number of simulations run, and the analytical computed probabil- 
ity of interaction 1 — exp(— x/A), where A is both, the one corresponding to 77 
absorption (left) and ICS (right). 

given the target and injection energy, the mean free path defines the distribu- 
tion for a sufficiently high number of runs. 

Basic geometry 

Parameters such as the viewing angle towards the observer, ao6s, the separation 
of the binary, and the distance to the shock region along the orbit, (Eq. 
[1]), are directly connected to 7-ray observational results and are thus crucial 
for detailed model discussion. Table H] gives the set of LS 5039 orbital and 
binary parameters that are relevant for modeling. These parameters come 
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Table 4 

LS 5039: system parameters 



Parameter 



Symbol Adopted value 



Radius of star 




9.3 i?0 


Mass of star 




23 Mq 


Temperature of star 




3.9 X 10"^ K 


Mass loss rate of star 


M 


10-'^ Mq yr^i 


Wind termination velocity 




2400 km s"^ 


Wind initial velocity 




4kms~^ 


Distance to the system 


D 


2.5 kpc 


Eccentricity of the orbit 


£ 


0.35 


Semimajor axis 


a 


0.15 AU ~ 3.5/?^ 


Longitude of periastron 


OJp 


226° 




Fig. 21. The geometry of the orbit of the inclination i in the binary system (with the 
MS - nia,ssiv6 stcir in the ccntGr) wliGrc the Sjiiglc to th.6 obscrvGr, o^ofos? is defined. 
No is the normal to the orbital plane, 9t is the angle related to the orbital phase 
and 0[ is its projection in the plane of the observer. The angle (3 gives the height of 
the pulsar for given phase above the observer plane. 

from the recent work by Casares et al. (2005b), either directly from their 
measurements or from the values compiled by them. The inferior conjunction 
(INFC) is defined as the phase when the binary is viewed from the pulsar side 
(the smallest aobs)- Exactly in the opposite side of the orbit we find superior 
conjunction (SUPC), when the pulsar is behind the massive star (the largest 
Oiohs)i see also Fig. [H One can see that for LS 5039, INFC (0 ~ 0.72) is close 
to apastron phase, while SUPC (0 ~ 0.06) is close to periastron. 



Figs. [21] and [22] shows the angle to the observer, function of orbital 

phase for two considered inclination of the binary, i = 30° and i = 60°. In 
addition, the separation of the binary is marked. The viewing angle changes 
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Fig. 22. The LS 5039 pulsar angle to the observer as a function of phase along 
the orbit for two different values of the binary inclination angle i. INFC, SUPC, 
periastron, and apastron phases are marked. 

within the limits (90° — i, 90° + i) and depends also on the longitude of the 
periastron Up. 

The comparison for the two inclinations shows already that aobs varies sig- 
nificantly in the case of the larger angle, as the difference is 6aobs = 120°, 
what is crucial for all angle-dependent processes discussed in this paper. The 
separation of the binary is given by 

- " (12) 



1 + e cos 9 



where e is the eccentricity of the orbit, 

the value of p is defined as 

p = a{l-e^) (14) 



and where b is the semi-minor axis. As the semi-major axis of LS 5039 is 
a = 0.15 AU = 3.5 Rs the assumed pulsar in the orbit is at periastron only 
2.25 Rs from the massive star, whereas at apastron, it lies at = 4.7 Rg, about 
a factor of 2 farther. 

To get the relation for the angle to the observer we use the geometry shown in 
Fig. [2T] where starting from the spherical triangle containing the inclination 
angle i, we have the relations: 

sin 9t cos i = cos /5 sin 9[ (15) 



38 



and 



cos 6t = cos P cos 9[ 



(16) 



which is the cosine theorem for a rectangular spherical triangle. Then, the 
angle to the observer defined in Fig. is calculated from the formula: 

cos^ P = sin^ 9t cos^ i + cos^ 9t, (17) 



from which: 

dobs = 11/2 ±13. 



The sign in this relation depends on the orientation to the observer, i.e., if the 
pulsar is below or above the orbital plane. Note that Eq. (fT7|) is not fulfilled 
all around the orbit as in here the angle 6t = denotes the common point for 
the orbital and observer plane. The relation can then be used after defining 
the orientation of the orbital plane with respect to the observer given by the 
inclination i and the periastron position ujper- The latter is done when we 
tilt the inclined orbital plane with Uper + 7t/2. To find the internal phase 6t 
dependent on the true anomaly 6 (with = at periastron), we have to find 
the spherical triangle from Fig.[2T]in the system of the orbital plane. Including 
the tilting of the systems we have the relation: 

et=e + UJper + TT/2. (19) 



Then we have to check if the calculated value of 9t is within one orbital phase 
< < 2n. If not the angle have to be replaced by 9t + 2n (if 9t < 0) or 
9 — 2n {9t > 27t). To use the Eq. (IT7I) and calculate the corresponding angle to 
the observer we have to recalculate the angle 9t once again as it is only valid 
in the range < 9t < 7i/2. We can divide the orbit in four quarters in which 
the following transformations are needed: 

(1) n/2 -9tii0<9t< n/2, 

(2) 9t - n/2 if 7r/2 <9t<n, 

(3) 37r/2 -9tiiTT <9t< 3n/2, 

(4) 9t - 3n/2 if 37r/2 < 9t < 2tt. 

Inserting such phase into the Eqs. (fT7|) and (fT8|) we get the angle to observer 
corresponding to the true anomaly 9. 
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Fig. 23. Basic parameters for photon and electron propagation and their interaction 
in the anisotropic radiation field of the massive star. The injection is produced at 
a distance Xi and angle a. Rg denotes the massive star radius, and the tangent 
direction to its surface is given by angle f3s- The propagation path is denoted by 
and Xe, while the interaction occur at distance xq. The scattering angle for 
electron-photon (e) or photon-photon (e) interactions is 9. 

Thermal radiation 



The source of the thermal radiation field is the massive star of early type (O, 
Be, WR). The spectrum is described by Planck's law, which differential energy 
spectrum (the number of photons of given energy e per unit energy de, per 
unit solid angle fl, per unit volume V) is given by: 

dn{e,Q) An 

"^'^ " dedQdV ^ {hcfe^/'^^^-r ^ ' 



where e is thermal photon energy, h is the Planck constant, and k is Boltzmann 
constant. 



Gamma ray absorption and e+e production: opacity and geometry 



The optical depth to 7-photon absorption in the radiation field of the massive 
star up to infinity can then be calculated from the integral: 

00 

T^^{Ej,Xi,a) = \~}^{Ey,Xi,a,Xy)dx^, (21) 



where is a photon interaction rate to e~^e~ production in an anisotropic 
radiation field and x^ is its propagation length. When the propagation occur 
toward the massive star surface the integration is performed up to the stellar 
surface. The photon interaction rate, A~3, is related to a photon of energy E.y 
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injected at a distance Xi from the massive star, at angle a (see Fig. [23]), and 
is given by the formula: 



''^{E^, Xi, a, Xy) = J (l + n)dn J ^(1) J n{e)a^^{P)de, (22) 



where is the distance to the interacting photon from the injection place 
along its propagating path. The integration limits will be discussed in a dif- 
ferent part of the Appendix. For now we focus on the variable /i in the first 
integration. It is the cosine of the photon-photon scattering angle fi = cos 6 
(see Fig. [23]) . The angle is the azimuthal angle between the photon (7-ray) 
propagation direction and the direction to the massive star, while 0^ gives 
its limit value for the direction tangent to the massive star surface (see Fig. 
^51) . The cross section to e'^e~ production is denoted {(3), where the 

parameter f3 in the center of mass system is 



(3 = Vuj^ - myu, (23) 



with 

1 



= -E^e{l + cose) (24) 



being the photon energy squared (in this notation it is assumed also that 
c = 1). From this we get 

= 1- 2m^/E^e{l + /i). (25) 



The kinematic condition for the angle 6 which defines the threshold for the 
e'^e~ creation is given by expression: 

2m2 ^ ^ 

> l^lim = -f^ 1- (26) 



To simplify the equations we rewrite the internal integration in Eq. fl22l) mak- 
ing use of hifi), such that, 

= [ n{e)a^^{P)de. (27) 



With the replacement 
^2 = 1- a/e. 
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where 



a = 2m^/E^{l + fx), 



(29) 



and defining the constant S = Sn/i^hc)^, the spectrum of thermal photons is 
now given by the formula: 



n(e) = S 



1 



(1 _ ^2)2 g(a/[(l-/32)fcT,)] _ l' 



(30) 



Substituting in the internal integral /i(yu), and introducing the integration 
variable to (3, via 



de = 2al3/{l- (3yd(3, 



(31) 



yields to the integral: 



/CL 1 
(1 _ ^2)2 e(a/{l-/32)fcT.) _ ^(^^^WP- 



(32) 



The lower limit of integration is from the energy condition for the process 
7 + 7 — >■ e~^e~ , i.e., it follows from the threshold condition E^e{l + fi) = 2m^, 
where uj = m. . The upper integration limit comes from the relativistic limit 
uj ^ m, where we get [3 ~ 1. To proceed forward, we introduce a dumb 
variable, 6, by a = /cT^fe, to get: 



h{h) = 2S J{kT,fb'a,,{(3)—^ 



P 



(1 - /52)4 e(V{i-/3^)) - 1 



dp. 



(33) 



The cross section for e~^e pair production is (Jauch and Rohrlich, 1980): 



^77(/3) = 2^0^(1 - r )[(3 - PI 1^ - 2/5(2 - n] 



(34) 



where tq is the classical electron radius, and ax = iTrrQ is the Thomson cross 



section. When putting this expression into the integral Ji (6) (Eq. [331) we get 
finally, 



(3-/3^)lni±|-2/3(2-/32) 



X 
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log(b) 

Fig. 24. The plot of the function C{b) = Ii{l3)/Ci, where the integral Ii is given by 
Eq. §5i). 



(3 



-d/3, 



(35) 



(1 - P^f e(V(i-/3^)) - 1 
with Ci = lQTi{kTslhcfaT. 

Weparametrizee the internal integral and write it as a function of 6, 

C{h)=h{P)/C^. (36) 



The plot for this function is shown in Fig. (12^ . Then, the integral we are after 
can be written as Ii{h) = CiC{h) and 

A-y^~^(-E'^, Xj, a, x^) = Ci / (1 + /i)c//i / d(j)C{h). (37) 



In the second internal integral, we have to fix the limits (having in mind that 
the angle depends on the angle 6, then also on the variable //), 

<t>s 

h= j C{h)d(l) = 20,C(6) = C(6)<l>(/i). (38) 

— Ips 

The angle 0^ determines the maximal azimuthal angle of 7-ray photon prop- 
agation with respect to the direction of the thermal photon, so that it gives 
the directions tangent to the star surface. This condition for 0^ can be deter- 
mined from the spherical triangle (shown in Fig. [251) . From the cosine theorem 
of spherical trigonometry 

cos (3s = cos 9 cos (vr — a) -|- sin 9 sin (vr — a) cos 0^, (39) 
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Fig. 25. Left: Spherical geometry for a photon absorption process in the radiation 
field of the star (a thermal photon is denoted with e). The angle Ps is defined by 
the line to the star center and tangent to its surface. The angles a and define the 
direction of gamma-photon incoming towards the star. The angle a is measured from 
the central line (joining the star and injection place) whereas the angle 6, from the 
line of gamma photon propagation. The directions of gamma-photon propagation 
defined as "inward" and "outward" with respect to the star are shaded. Right: two 
dimensional geometry for definition of the parameter // = cos 9, where ni = cos 9f 
and iJ,2 = cos 6*1 are given by Eq. ([^2]) . The azimuthal angle (ps is the angle between 
the direction of gamma propagation and propagation of e.. 

where Ps is defined as in Fig. (125|) . After some algebra we get: 
cos /5s + /i cos a 

COS0, = . , (40) 

V 1 — sm a 



where sin (3s = Rs/^i, while < 0s < ^- The angle a determines the direction 
of incoming 7-ray photon (see Fig. [2511 . 

The limits of integration with respect to the parameter /i (see Fig. [23]) . are 
the range of angles for soft photons coming from the star, 

X^^~-^{E^,Xi,a,x^) = Ci J {1 + dfiC{b), (41) 

Ml 



where /ii = cos 61, /i2 = cos 6*3. Using basic trigonometrical dependencies be- 
tween the angles, Of = n — a — P and 62 = 7!' — a + P (Fig. [25], left), we get 
two expressions: 



Hi = — cos a cos Ps + sin a sin Ps 
fj,2 = — cos a cos Ps — sin a sin Ps 



(42) 



Depending on the angle of 7-photon propagation a we can precise the condi- 
tions for the above integration limits. We distinguish three regimes: 

(1) the range for a towards the star surface: {n — Ps < a < n), 
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(2) the "middle" range: {Ps < a < n - Ps), 

(3) the range outwards of the massive star surface: (0 < a < f3s)- 

Integration over /i can be simphfied by separating the main integral, Eq. (I4ip . 
in the sum of the corresponding integrals for these regimes, each depending 
on their limits fJ^i, fi2, fJ'iim (see Eq. |12] and [26]). A simplification occurs for 
the specific range of the angle 6 for which the function = 2(j)s (Eq. 

I38l) is maximal, i.e., $(/i) = 27r, what corresponds to the scenario when the 
photon propagate toward the star. Another simplification happens when $(/i) 
is minimal, i.e., $(/i) = 0, for an outgoing photon. To simplify the formulae for 
numerical implementation we introduce the function f{fi,b) = CiC(6)(l + /i) 
such that 

X-y^~^{E^,Xi,a,Xy) = J f{jj,b)^{iJ?) dfi, (43) 

Ml 

where we introduced internal function f{fi, h) = Ci(l + /x) C{b). To fasten the 
simulations the integration limits can be fixed following the rules (note that 
relation (/xi > fj,2) is always fulfilled): 

(1) for "inward" propagation : 

• if {fJ'lim < /i2 < /ii): 

^77~^(^7' ^7) = /A' /(/^' ^)'^(/^) ^/^ + /(/^' ^)27r c^/i, 

• if (/U2 < fJ'lim < fJ'l)- 

• if (yU2 < yUl < /iiim): 
A^^"^(E^,Xi,a,a;^) = J^^^^ f{fi,b)2ndfi. 

Apart from these, the case of a = vr can also be separately defined (incom- 
ing at a central spot of the star) when {fii = 1^2) and \^^~^{E^, xi, a, x^) = 
IL^n /(/^' b)2ndfx, where the smaller value from {fii,fiiim)- 

(2) for the "middle" range of propagation: 

• if if^lim < f^2< fJ'l)- 

\7~^(^7' ^7) = /^a' /(/^' ^)'^(/^) ^/^' 

• if {f^2 < fJ'hm < /Ui): 

X^^'\E^, Xi, a, x^) = Jl^^^^ fifi, 6)$(/i) dfi, 

• if (/i2 < yUi < /iiim): 

(3) for "outwards" directions: 

• if {fJ'lim < /i2 < /ii): 

^77"^ (^7' 3;^, a, Xj) = Jl^^ /(/i, &)$(/i) c//i + J;;^'^ /(/i, 6)27r c?/i, 

• if (/i2 < IJ'lim < /ii): 

\7~^(^7' 3;^) = /wL /(/^' ^)'^(/^) 

• if {1^2 < /^i < yWiim): 
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Apart from these, as before, the case of a = can also be defined (es- 
caping radially from a central spot of the star) when {ni = 1x2) i if 
{ill > fxiim), and we have X^^~^{E^, Xi, a, x-y) = /^|^ /(/i, 6)27r d/i, other- 
wise Xyy~'^ {Ej, Xi, a, x^) = 0. 



The energy spectrum for e^e pairs produced by 77 absorption process 



To simulate the energy of an electron (positron) produced in the process of 
gamma absorption we follow the energy distribution of e~^e~ pairs (see Eq. 
[TTl) . The e~^e~ pair spectrum produced by a photon of energy at a specific 
distance Xp is given by the expression: 



dW{E^, Xi, a, Xy 
dEedxy 



y (1 + n)dn j d(t) j n{e)a{Ey, E^, n)de, (44) 



where E^ is the energy of the interacting photon, x^ is the propagation length, 
Xi the distance of the injection place form the massive star, a is the angle of 
propagation, and E^. is the energy of produced lepton (e+e~). The limits of 
integration are discussed in the following paragraph. 

The energy spectrum of thermal photons e is again given by the Planck's 
spectrum n{e) (Eq. [20]). The cross section for pair production a{E^, Eg, fi) 
depends on the photon energy E^ and energy of electron E^. The variable 
/i = cos^^, angles 9 and 0^ are defined in Fig. [25ll . 

The cross section is given by the formula (Akhieser & Berestezki, 1965): 



^ " ^' 2 E, e,/5,7c ^ m^B\E,, e,) + elA\E,, e,) 
2{el-m'fA\E,,e,) 
{m^B\E,,e,) + elA\E,,e,)Y^' 



(45) 



where Ec is the photon energy in centrum of momentum (CM) of the interacting 
photons el = ^Eje{l + /i). The parameter 7c = (-E^ + e)/2ec is Lorentz factor 
in CM system and m = meC^ is the electron mass. In Eq. fHSj) . the following 
functions are introduced: 



{Ee - -fcSc 



(7c/5c)2(£2-m2)' 
A'{Ee,ec) = l-B\E„e,), 



(46) 
(47) 
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Fig. 26. The integral function D(X) n) (Eq. [60]l . where x = 2kTsEy[l + fi)/m'^ and 
K = Ee/E^ are the parameters. 



where f3c is the velocity of the CM system in units of the speed of hght, 



7c/3c ~ \Jlc ~ "^^^ internal integration in Eq. fl44p . i.e., that performed over 
the thermal photons spectrum, according to Eq. fHSl) . is: 



in 



{hey J {e^/'^Ts _ 1) e^p^^^ 

m^B^{E,,e,)+elA^{E,,e,) 
2{el-m')'A\Ee,e,) 



{m^B\E,,e,)+elA\E,,e,)Y 



de. 



(4J 



The lower integration limit is given by the condition: 



m?E^ 



2Ee{E^ - Ee){l + 12) k{1-k) 



(49) 



In the last expression, we have introduced the new variables n = E^/E^ and 
F = w?/2E^{l + /i), where k E (0, 1). 
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For a photon energy ^ e, we have 7c ^ 1, — 1, and ScPcIc = 
After rewriting Eq. (HHj) we get 

8 2 °° 

= / (50) 



where, to simphfy the formula, we have defined 



fie) 



m^B^{E,,e,)+elA\E,,e,) 
2{el-mYA\E,,e,) 
{m^B\E,,e,)+elA\E,,e,)y 



(51) 



To proceed forward, we introduce the variable z = e/kTg, so that kTgdz = de, 
and the integral is now given by the expression: 



oo 



(52) 



where the lower limit of the integral is given by Zmin = ^min/kTs, 

Zmin = m^/{2E^kTs{l + h)k{1 - k)) = m^/xK{l - k) (53) 

and a new parameter, 

X = 2E^kT,{l+fi), (54) 

was introduced. The parameter x is limited by the condition > rn^, so then 
X > ^rr? jzmax = 4:kTsm^/emax, (55) 



where emax is the maximal thermal photon energy from the blackbody spec- 
trum (we assume emax = SOkTg). By using the variable 

g = el = ^E^e{l + /i) = zx/^, (56) 



the additional needed functions are given by expressions: 
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2. ^ g{2K-iy 
B (x,«:) = 5—' (57) 

A\x.t^) = l-B\x,t^)- (58) 

Function (!5T|) can now be rewritten as: 



^^^^ ^ 2g-m^ + (g-m^)A^(x, /t) 
2ig-myA\x,f^) 



2- (59) 
m^B^{x, /«) + 0A^{x, /«)] 



We can further conveniently replace the parameter x with a dimensionless 
form X = 2kTsEj{l + 11) /w? to get finally the full set of needed magnitudes 
in a useful way for integrating: 

h{E,,E,,^i) = —±-^—D{x,k), (60) 

00 

D{x,^)= J ^^f{z)dz, (61) 

^ 2g-l + (g-l)A2(x,/t) _ 2{Q-lfA\x.^^) .g2) 

^^""^ 52(x,/€) + ^?A2(x,/€) [52(x,K) + f?A2(x,K)]'' ^ ^ 

A\x,'^) = 1-B'{x,k). (64) 



To produce the electron of the energy k the process is limited by the condition 
z > Zk, where 

^ (65) 



x/«(i - 



On the other hand, from the kinematics of the process we have limitation 
z > z^, where 

_ 4 _ 2m^ 



The lower limit of the sought integral is the greater value of the two expressions 
quoted above. The integral given in Eq. (!6Ti) can be tabulated with respect to 
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the parameters x ^i-nd k. The plot of the function D{xi i^) (Eq. [6T]) is presented 
in Fig. ([26D. 

The remaining internal integral in Eq. (jUj) is given by the function = 
/-f^^ dcf) = 20s, where the angle 0^ is defined by Eq. ( HOl) . With this in mind, 
the formula for the spectrum of pairs is given by: 



dW{E^,Xi,a,x.y) C2 



J{l + fiMfi)D{x,K)dfi, (67) 



where the constant C2 = ScT'rTn'^ (kTs)'^ / {hcY . The integral limits, {fii, (12) can 
be calculated according to Eq. following the rules described in previous 
paragraphs. 

The e+e^ pair spectra as a function of initial parameters for the process. 
Fig. I23p . and the energy of interacting photon, E^, are shown in 
Fig. fl27p . The produced e+e^ pair energy spectra are symmetric with respect 
to the electron energy E^/E^ = 0.5, what is characteristic for the gamma 
absorption process. Figures show also that the process of pair production 
strongly depends on the assumed geometry of interacting particles so that 
the efficiency is decreasing for gamma photons propagating in the outward 
direction and if injected farther from the star surface. Note that the presented 
spectra are calculated for specific parameters of the massive star, but serve as 
a general example. 



The production of high- energy photons in Inverse Compton scattering 



The electron propagates in the radiation field of density n{e), and the electron 
injection place is defined by the distance Xi and angle a with respect to the 
massive star (see Fig. [2^ . The thermal photons are again described by a 
blackbody energy spectrum (Eq. [20]) . 

The rate of interaction of a high-energy electron of energy E^, at the distance 
Xe from its injection place, is given by the integral: 

\-Wz7 ^ ^ } dN{Ee,Xi,a,Xe) , . 

Xjcs[Ee, Xi, a, Xe) = - J d^^t ^ ' 





The integration of the scattered photon spectra dN{Ee, Xi, a, Xe)/dE^dt is 
taken over the energies of the produced photons, E^, and is limited by the 
maximal photon energy scattered E™'"'^. 
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E^=10^MeV 




Fig. 27. The differential spectra for e^e~ pairs produced in the absorption pro- 
cess of a photon with energy E-y. The massive star has here a surface tempera- 
ture Tg = 10^ K, and radius Rg = 10 i?©. The spectra for a fixed place of pho- 
ton injection, defined by Xi = 2Rs, and a = 170° for different energies of pho- 
ton are shown at the left in the top panels. The energies of gamma photon are 
= 10^, 10^ 10^, 10^ MeV, from top to bottom curves. The same dependence but 
for Xi = lORs is shown to the right in the top panels. The e~^e~ spectra from 
different angles a at the same injection distance Xi = 2Rs and energy of photon 
E-y = 10^ MeV are presented in the left, bottom panels. The spectra are calculated 
for a = 30°, 60°, 90°, 120°, 150°. The dependency on the distance of the photon in- 
jection place, for Xi = 2, 5, 10, 15, 20, 30 Rs , and fixed E^ = 10^ MeV and a = 170° 
are shown in figure to the right, bottom panels. The energy of leptons are normalized 
to the injected photon energy. 




Fig. 28. The geometry for the Inverse Compton scattering in the observer (LAB) 
frame. The Lorentz factor of an electron is 7e, the energy of a photon before scat- 
tering is ai, and angle of scattering is 6i (the angle between the directions of prop- 
agation of the low energy photon and the relativistic electron) . 
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Fig. 29. The kinematics of IC scattering in the center of mass system (CMS). The 
photon before scattering is denoted by a[ = e/nieC^, and after it, by a'2 = E^/nieC^. 
The propagation of the photon before and after interaction are given by angles 6[ 
and 02, both are polar angles measured with respect to the velocity vector of the 
electron /3 = v/c. The scattering angles are the polar angle x\ the azimuthal 
angle 4>'. 

We calculate the optical depth to IC scattering, at the propagation path x^, 
integrating the rate of electron interactions in the anisotropic radiation field: 



If the electron propagates directly towards the massive star, the path is ter- 
minated by the massive star surface. 

The kinematics of Inverse Compton scattering 

We adapt here the formulas from Jones (1968) (therein Eq. 1 through 10) and 
Blumenthal and Gould (1970) (their chapters 2.6-2.7), wherein the scattered 
photon spectra are calculated in the isotropic radiation case. In our work, 
we integrate over the angles of incoming low energy photons (the solid angle 
which defines the surface of the massive star) what differs from integration 
over whole solid angle in the isotropic radiation field. 

The energy of photon before scattering in the LAB system (see Fig. [28l) is 
denoted by ai = e/mec"^, while that after scattering, by 02 = E^/nieC^. In 
the CMS (see Fig. f2U\i all quantities are marked with prime variables, so that 
the energy of the photon before and after scattering are written as a[ and 
The angle between the direction of the electron propagation and the incoming 
stellar photon is denoted as di {6[ in the CMS), while the angle between the 
electron direction and the photon after scattering, is 62 (^2 CMS). These 

are polar angles measured with respect to the electron velocity vector = v/c. 
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The angles of scattering in the electron rest frame generate a spherical triangle, 
so that 



^2 = cos 9[ cos x' + sin 9[ sin x' cos 



cos 62 = cos 9[ cos x' + sin 9[ sin x' cos 0' (70) 



(see Fig. [29]). The energy of the scattered photon in the CMS is given by the 
formula: 



l + aUl-cosxO' ^^^^ 



where x' is the photon scattering angle, defined as the azimuthal angle between 
the propagation direction of the photon after scattering, and the photon 
propagation direction before it, a[. 

The energy of photon in the LAB system can be found from a Lorentz transfor- 
mation. In a simplified case, when the photon propagate head-on with respect 
to the electron direction, and after interaction, it has the energy: 

"2 = 7ea2(l + /5cos(7r - x')) ~ 7ea2(l - cosx')- (72) 

In the electron's rest frame (see Fig. 1291) . for the polar angle 6[ = we have 
x' = (see Eq. [72]), what means that the photon energy after scattering (see 
Eq. [72]) does not depend on the azimuthal scattering angle (J)'. 

We can calculate the maximal energy of the scattered photon from relations 
( ITT]) and ( !72l) . such that amax = 27ea2- ^^e Thomson limit {a[ <^ 1), the 
maximal photon energy depends on the electron energy and equals a^ax ~ 
47^01. In the Klein-Nishina limit, which is the highly relativistic case, this 
energy is proportional to the initial electron energy and it is a^fl ~ 7e«i- 



The IC scattered photon spectrum in an anisotropic target field 



In the case of mono-energetic isotropic photon background (in the LAB sys- 
tem) the angular spectrum of incoming low energy photons in the electron 
rest frame is given by the formula (Jones 1968): 



n {e^)dcose^ = — — — 

27e(l — p COSb'^ 



(73) 



For relativistic electrons (7e ^ 1), the incoming target photon polar angles 
are in the range < ^^'^^ < 6*^^25 where 6*1/2 ~ l/7e is the angle of a photon 
'cone' apex. 
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0.25 



0.5 



0.75 



Fig. 30. The photon spectra (Eq. [951) for fixed initial parameters (the electron in- 
jection place): the distance to the massive star Xi = lORg and angle a = 180° (see 
Fig. (I23p ). as a function of electron energy E^. The temperature of the star is herein 
taken as = 10^ K and the radius Rs = 10 Rq. The electron Lorentz factors shown 
are : je = 10^ 10^, 10^ 10^ 10^ 10^ 

The photon energy spectrum before scattering in the electron rest frame is: 



for 7e ^ 1, where function Q{x;a,b) is the characteristic function defined in 
a range (a, b) and takes the values: 



The cross section to IC scattering in the Klein-Nishina range is given by the 
formula: 



where y' = cosx', ''"o = e'^/mc^ is the classical electron radius and the internal 
function is given by 




(74) 





5(4 -/(«;, y')) 



(75) 




(76) 
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The number of scatterings per unit time t' is n'ca. If dN/dt — ^dN/dt', 
then after integration over (j)' one gets: 



TTTqC 



l+y 



n 



dtda[da'2dy' 2ah^ [1 + a[{l - y')]^ 



X 



1 + 



X 



From the relation 

da'^da^dy' = [1 + a'^{l - y')f da'^dy' df 



(77) 



(78) 



and after integration over / we get: 



dtdai^dy' '^Oilj^ 



.l\2 



e L 



l-a'^il-y') 



X 



X 



l-«2(l-?/0 



(79) 



Applying the Doppler-shift formula to the energy and LAB frame energy 
^2 yields to relation 



^2 = Q;2/7e(l - /?y')- 



(80) 



Then, by means of the replacement p—1 — Py': 



d^N 



dtda2dp 2a\^^{l — Q;2/7e) 



7rroCQ;2 



X 



1 - 0:2/76 



6(p; Pi,P2) 



(81) 



which is obtained under the assumption that (1 — y') ~ p. The limiting values 
in the argument of the function ©(p; pi, P2) are: 



02 



P2^ 



20172 (1 - Q;2/7e) ' 
2^2 

ai(l - a2/7e)' 



(82) 
(83) 
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where the parameter p is hmited by conditions pi < p < Pp, with 



Pi 



Pp = 2 



To calculate the spectra of scattered photons from the spectrum (Eq. [8T]) as a 
function of our parameters, dN/dtdE^{Ee,Xi,a,Xe), we replace = a2'meC^ 
and e = aime(? and integrate over the background photon spectra n{e) and 
parameter p, 



dN 
dtdE. 



27f (mc2)3(l - E^/Ee) 



$(/"(p)) 

Q(p;pi)P2) 
p2 



p2 - 2p + 2 + 



dp, 



^(ei) , 
— T^dei X 

1 — E^J Ee_ 



X 



^5) 



where 



p = E^mc^ / {e'je{Ee - E^)p) - 1, 



^6) 



and the function $(p) has been already defined in previous paragraphs as 
For convenience, let us rewrite the internal integral on p as follows: 



Pmax 

G(E^,e,p)= J $(p) 

Pmin 

0(p;Pi,P2) , 

2 

p2 



p^ - 2p + 2 



1 — E^/ Ee 



X 



^7) 



If the range of parameter for p from Eq. (IHTI) is consistent with the range 
{pu Pp) (Eq. [HI]), then the integration limits {pmin, Pmax) correspond to the 
common part: {pi, P2) H {pi,pp), and the integral over p is given by function: 



G{E^, e, p) = 27r 



2 {E, E,)^ 

p-21ogp ^ 

p (1 - E^/E^)p 



pmax 
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/ E, / E, 

Fig. 31. Left: The photon spectra (Eq. ([96|) ) for a specific electron energy 7e = 10^ 
and fixed distance to the massive star Xi = WRs {xe = 0), calculated for different 
electron propagation angles: a = 0°, 30°, 60°, 90°, 120°, 150°, 180°. Right: The pho- 
ton spectra depending on the electron energy Eg and the initial injection parameters. 
The spectra are calculated for 7e = 10^, the angle of propagation a = 180° and sub- 
sequent distances to the massive star: Xi = 5, 10, 20, 30, 50 Rs (vs = Xi/Rs,Xe = 0) 
(right panel). The temperature of the star is herein taken as Ts = K and its 
radius is Rg = 10 Rq. 

From the Doppler formula and transformations between the systems we get, 
p = E^/72e(l + /3cos^2)(l - E^/E,), (89) 



so that the integration limit over p, for angles Of and 62 (see Fig. [25!) are given 
by the formulas: 



Pi(^i) = -7r-^T7Z^h^^rT^^ (90) 



e7e(l + /?cos^f)(l- 
Ej 

e7e(l+/?COS^i)(l- V^e; 



Pi02)= , . ^ (91) 



where the angles 6f, 62 are the angles which defined the directions of low 
energy photons coming from the star. In addition, the parameter p is limited 
by kinematic condition for the scattering to happen, 

cosekin> 2ml/ E^e-1. (92) 



The integration in respect to p can be separated to get the sum of integrals, as 
for p > p{9kin) the internal function $(/u) = 27r, and the function G{E^, e, p) 
is given by (IHHj) . If the function $(/i) = 2tt we can integrate the function 
G{E^,e,p) (Eq. [87]). Otherwise, if $(/i) = 20^ we can rewrite the function 
G{E^,e,p) in a compact way as follow. Thus, these cases are summarized by 
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ymax 

G{E^,e,ii)= j 2^sGp{E^,e,ii)dp, (93) 



Pmin 



G{E„ e, /.) = 271 [GUE„ e, . (94) 



where G,{E^,e,fi) = [p" - 2p + 2 + {E^/ E,f / {I - EJ E,)] x e(p; pi, pa)/^' 
and Gint{E^,e,p) = p-2\og p-2/ p-{E^ / E^f / {{I - E^/ Ee)p). The hmits of 
integration depend on combinations of {p{6i), p{62)) and {pi,pp), and respect- 
ing the constraints put by p{6kin)- The final integration G{E^, e, p) is given by 
the formulas: 

(1) for "outwards" propagation : 

• if {Pmax < Pi): G{E^,€,p) = 0, 

max 

> pi) and {pmin < Pi) ■ 
and {pmax < P2) ■ 

G{E„ e, p) = 2^,Gp{E^, e, p) dp, 
or {p 

max > P2) ■ 

G{E^, e, p) = 20,G',(E^, e, p) dp + 27r [G,nt{E,, e, /x)];^"% 

• if {P max 

> pi) and {pmin > Pi)' 

if {Pmin < P2) and (p^a^: < P2): 

G(E^, e, p) = 2^,Gp{E^, 6, p) dp, 

if (Pmin < P2) and {pmax > P2) ■ 

_ G(E^, 6, p) = /;^^^ 20,G,(i?„ 6, p) rfp + 27r [G,,,(i?„ e, p)Y-% 

if {Pmin > P2)' 

G{E^, e, p) = 27r [Gi„t(£;^, e, p)]^"^ 

(2) for the "middle" range of propagation: 

• if (P max 

< Pi) or {p^in > P2): 
G'(E^,e,p) = 0, 

• if (P max 

> pi) and (p 

max < P2): 

G(i?^, e, p) = /rAl(p.,„,pO 20.G',(£;^, e, p) dp, 

• if {pmax > Pi) and (p 

max > P2): 

G{E^, e, p) = I^AX(p^,„,p^) '2(psG^{E^^, e, p) rfp. 
where MAX [a, b) gives the larger number from the brackets. 

(3) for "inward" directions: 

• if {pmin > P2) ■ G{E^,e,p) = 0, 

• if (P min 

< P2) and {pmin > Pi)' 

if (P max > P2) ■ 

G{E„ e, p) = 2<P,G,{E^, e, p) dp, 

if (P max < P2) '■ 

G{E^, e, p) = 2(t)sGp{E^, e, p) dp, 

• if {Pmin 

< P2) and (p^in < Pi): 
if (P max > P2) : 



G{E^, e, p) = /;f 2<PsGp{E^, e, p) rfp + 2n [Gint{E^. e, p)] 



pi 



if (P max 

< P2) and (p 

max > Pi) : 
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if {Pmax < P2) and (p 

max < Pi) ■ 

Finally the rate of electron scattering on its propagation way to IC process is 
given by expression: 



^max 



dN 
dE^dt 



{Ee, Xi, a, Xe) dE. 
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(95) 



where, 



-{Ee, Xi, a, Xe) = C ^ 



dtdE^^ " " ' -f^{l-E^/Ee) 

^^j^^G{E„e,p), (96) 



with the constant C = rQc/2(mc^)^. The function G{E^,e,p) is described by 
Eqs. ( IHTl) and (IHHl) . which have to be chosen with respect to the value of the 
function $(yu) and the limits to the parameter p, as explained. The lower limit 
for the integral in Eq. fl96|) corresponds to the lowest energy of the scattered 
photon 

emin = E^/4jI (97) 



while the upper one in the integral of Eq. (!95l) is the maximal scattered photon 
energy 

„ = 4^2^. (98) 



The photon spectra calculated from Eq. fl96|) are presented in Figs. fl30l) and 
( I3TI) . They are calculated with respect to the electron energy Ee and its initial 
parameters of injection close to the massive star (see Fig. [23]). The photon 
energy spectra reflect the features of the cross section for IC scattering, with 
characteristic peaks at the high energy range (for E^/Ee ~ 1) - Klein-Nishina 
regime if the energy of the incoming electron is relativistic. These plots shows 
also the dependence on the propagation angle (with the highest production 
rate for direction toward the massive star) and place of the electron injection. 
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